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INTRODUCTION 


This  report  describes  theoretical  work  aimed  at  estimating  the  dis¬ 
persion  of  oil  engendered  by  the  joint  effects  of  turbulence  and  the  large 
scale,  wind-generated  convective  motions  in  the  ocean  known  as  Langmuir 
circulations. 

The  physical  problems  of  interest  in  oil  spills  concern  the  location 
of  spilled  oil  as  a  function  of  time.  Oil,  almost  always  being  lighter 
than  water,  tends  to  float  on  the  surface.  The  natural  action  of  water 
turbulence,  particularly  that  generated  by  breaking  waves,  is  known  to  be 
capable  of  disintegrating  floating  oil  layers  into  small  droplets,  which 
can  then  be  entrained  into  the  water  column.  This  process  is  dramatically 
accelerated  by  the  addition  of  chemical  dispersants  (McCarthy  eii  al, 
1978),  and  is  the  basis  for  a  method  of  oil  spill  cleanup. 

Although  most  attempts  to  model  the  distribution  and  disposition  of 
oil  in  the  sea  have  dealt  only  with  surface  oil,  vertical  dispersion  of 
oil  also  occurs  and  is  important  to  understand.  Since  mechanical  methods 
of  oil  spill  cleanup  rely  upon  oil  being  on  or  very  close  to  the  water 
surface,  the  ability  to  estimate  the  fraction  of  oil  driven  down  into  the 
water  column  as  a  function  of  sea  state  may  be  a  key  question  for  oil 
spill  clean-up  operations.  Furthermore,  oil  dispersed  by  chemical  means 
continues  to  be  buoyant  (Canevari  1978)  and  will  rise  to  the  surface  in 
quiescent  water.  Thus,  the  problem  of  vertical  dispersion  allows  one  to 
estimate  the  success  expected  when  dispersants  are  employed.  One  of  the 
advantages  of  vertical  dispersion,  whether  by  natural  processes  or  promot¬ 
ed  by  chemical  means,  is  the  division  of  oil  into  volumes  more  accessible 
to  biodegradation  (Canevari  1978);  the  vertical  extent  of  the  dispersion 
may  play  a  role  in  determining  which,  if  any,  biological  forms  are  able  to 
attack  the  oil. 

The  problem  of  the  vertical  dispersion  of  oil  has  not  been  exten¬ 
sively  studied.  The  first  published  work  seems  to  be  that  by  Leibovich 
(1975),  followed  by  Raj  (1978)  and  by  Milgram  et  al  (1978).  At  the  time 
those  attempts  were  written,  there  was  little  observational  evidence  of 
vertical  dispersion  of  oil  except  under  circumstances  (Forrester,  1971; 
Hess,  1978)  in  which  the  spill  occurred  rather  near  the  shore.  (In  the 
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latter  cases,  onshore  currents  and  surf  may  result  in  a  sediment  and  oil 
agglomeration  with  a  composite  density  greater  than  seawater.)  The  blow¬ 
out  of  the  Ixtoc  I  oil  well  in  the  Gulf  of  Mexico,  beginning  in  June,  1979 
and  continuing  for  more  than  four  months,  allowed  prolonged  observation  of 
an  oil  spill  clearly  away  from  surf  zones  and  hundreds  of  miles  from  its 
source.  Visual  observations  (made  by  divers)  of  subsurface  oil,  though 
few  in  number,  were  remarkable  (Williams  (1979),  Robinson,  Galt  (1980), 

Hooper  (1981)).  Apparently  oil,  in  the  form  of  flakes  the  size  of  corn 
flakes  and  thought  to  be  the  fragmented  semisolid  skin  of  large  oil  "pan¬ 
cakes",  was  mixed  to  a  depth  of  40  feet.  The  explanation  for  the  exist¬ 
ence  of  buoyant  particles  in  substantial  numbers,  originating  at  the  sur¬ 
face  and  mixed  to  depths  such  as  these,  defies  all  self-consistent  models 
of  turbulent  transport  previously  known  to  us.  It  seems  almost  certain 
that  the  explanation  must  reside  in  the  existence  of  a  mixing  mechanism  of 
large  scale. 

The  likely  mixing  mechanism  is  related  to  a  phenomenon  known  to  the 
scientific  world  for  more  than  a  century  (see  the  discussion  of  observa¬ 
tions  by  James  Thomson  and  his  brother  William,  Lord  Kelvin  (1862));  these 
are  windrows,  long  bands  of  flotsam,  or  compressed  organic  film  or  foam 
from  breaking  waves,  that  are  approximately  parallel  to  the  wind  direc-  ^ 

tic.  Windrows  are  now  known  to  be  the  visible  surface  manifestations  of 
large  scale  subsurface  vortical  motions  discovered  by  Langmuir  (1938),  and 
known  as  Langmuir  circulations.  The  rows  mark  lines  of  surface  conver¬ 
gence,  below  which  water  sinks  at  speeds  orders  of  magnitude  larger  than 
turbulent  f’uctuations  and  comparable  to  the  surface  wind  induced  water 
drift.  Since  oil  provides  an  ideal  visual  marker,  windrows  marked  by  long 
bands  of  collected  oil,  often  in  the  form  of  ropy  water-in-oil  emulsion 
referred  to  as  "chocolate  mousse",  are  much  more  clearly  visible  when  oil 
is  spilled  than  under  natural  circumstances.  They  have,  consequently, 
often  figured  in  accounts  of  oil  spills  -  the  earliest  one  we  know  of  is 
one  reproduced  in  N.K.  Adams  (1936)  report  to  the  Royal  Society  of  London 
on  oil-spill  problems.  More  recent  accounts  of  oil  bands  in  windrows  in¬ 
clude  Battelle's  (1969)  review  of  the  Santa  Barbara  blowout  and  Galt 
(1978)  in  his  discussion  of  the  Amoco  Cadiz  and  Hawaiian  Patriot  spills. 
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Accounts  of  the  windrows  in  the  Ixtoc-I  spill  indicate  that  they  were  a 
dominant  feature,  and  are  carefully  described  by  Atwood  et  al  (1980). 

The  simultaneous  appearances  of  a  banded  surface  structure  and  the 
appearance  of  oil  at  surprisingly  great  depths  are  quite  likely  related; 
the  working  hypothesis  of  the  investigation  reported  here  is  that  both 
features  result  from  motions  in  Langmuir  circulations.  If  so,  when  oil  in 
windrows  is  evident,  and  when  the  oil  present  can  be  broken  into  small 
particles,  then  it  is  likely  that  oil  particles  will  be  carried  down  into 
the  water  column  by  this  mechanism.  Our  object  here  is  to  try  to  estimate 
the  amount  of  oil  held  in  suspension  in  this  way,  given  the  availability 
of  standard  environmental  parameters  such  as  wind  speed  and  sea  state,  and 
perhaps  visual  observations  of  surface  windrows. 

The  complicated  processes  by  which  a  floating  oil  mass  can  he  broken 
into  small  droplets  or  fragments  is  considered  at  length  by  Milgram  et  al 
(1978),  and  we  do  not  address  the  question  here. 

While  the  deep  penetration  of  oil  is  attributed  to  persistent  large 
scale  motions,  the  possibility  of  any  oil  being  dispersed  into  the  water 
column  is  due  to  turbulence.  As  in  most  problems  in  fluid  dynamics,  the 
treatment  of  the  turbulence  is  the  most  difficult  and  uncertain  part  of 
the  problem.  A  purely  theoretical  approach  may  be  disasterous;  for  exam¬ 
ple,  reasonable  people  can  easily  arrive  at  estimates  for  the  dissipation 
rate  of  turbulent  energy  (a  key  figure)  that  differ  by  ten  orders  of  mag¬ 
nitude.  To  approach  the  question  of  dispersion  in  a  rational  way,  we  must 
introduce  the  best  available  information  about  turbulence  in  the  surface 
layers  of  the  sea.  Although  experimental  information  is  accumulating  at 
an  accelerating  rate,  data  on  oceanic  turbulence  is  not  plentiful  and  much 
of  the  existing  data  is  of  dubious  quality.  Fortunately,  pertinent  new 
turbulence  measurements  below  surface  waves  in  the  ocean  have  been  made  in 
the  last  two  years,  and  even  more  recently  in  Lake  Ontario,  which  provide 
important  guidelines  for  our  investigations. 

Assuming  that  turbulence  data  are  available,  it  remains  to  determine 
their  effect  on  the  transport  of  oil.  The  question  posed  in  the  present 
work  is  the  determination  of  the  dispersion  of  oil  due  to  mean  motions 
generated  by  Langmuir  circulation.?,  which  are  calculated  using  a  theory 


due  to  Craik  and  Leibovich  (1976)  (see  Leibovich  1980  for  a  summary),  and 
superimposed  turbulent  fluctuations.  A  model  for  the  motion  of  the  oil 
incorporating  turbulent  fluctuations  is  based  on  the  work  of  Lumley  (1978) 
and  of  Csanady  (1963)  and  is  not  straightforward;  and  a  considerable 
amount  of  attention  is  devoted  to  it  in  this  report. 

We  have  explored  the  joint  effects  of  Langmuir  circulations  and  tur¬ 
bulence  by  computing  the  distribution  of  oil  in  27  cases  illustrating  a 
range  of  surface  wind  speeds,  and  sea  states  and  for  several  oil  droplet 
terminal  velocities.  The  results  of  our  study  may  be  summarized: 

1)  Under  typical  sea  states  and  wind  conditions,  Langmuir  circula¬ 
tions  gather  oil  into  windrows  in  volumes  that  can  be  estimated 
by  the  method  of  §8. 

2)  Oil  may  be  trapped  below  the  surface  in  extensive  zones,  which 
we  call  "Stommel  retention  zones"  (or  "SRZ")  in  recognition  of 
earlier  work  by  Stommel  (1949)  on  a  related  problem.  Each  of 
these  zones  lies  beneath  a  windrow  and  its  size  and  placement  is 
determined  by  the  oil  terminal  velocity  and  the  Langmuir  circu¬ 
lation  Speeds.  Oil  particles  can  be  found,  with  equal  probabil¬ 
ity,  anywhere  within  the  zone.  The  size  of  these  zones  of  trap¬ 
ped  oil  easily  accounts  for  the  deep  penetration  observed  in  the 
Ixtoc-I  spill.  For  example,  for  one  illustrative  case  evaluat¬ 
ed  here  (oil  particle  terminal  velocity,  =  OJj  cm/sec.,  wind 
speed  =  15  m/s,  significant  wave  height  =  70.4  cm),  the  SKZ 
penetrates  to  a  depth  of  24  meters. 

3)  A  Stommel  retention  zone  is  fed  by  oil  transferred  from  the  win¬ 
drow  lying  above  it.  The  fraction  of  oil  trapped  in  a  SRZ  can 

be  calculated  by  the  method  of  §6.  The  fraction  of  oil  in  an 

SRZ  depends  upon  the  efficiency  of  transfer  from  the  windrow, 

and  this  depends  strongly  upon  the  distance  between  the  windrow 
and  the  top  of  the  SRZ.  This  distance,  in  turn,  depends  upon 
the  speed  and  length  scales  in  the  Langmuir  circulations,  and 

upon  the  terminal  velocity  of  the  oil  particle.  In  the  cases  we 
have  computed,  the  fractions  of  oil  trapped  exceeds  90%  for  ter¬ 
minal  velocities  as  large  as  1  cm  for  wind  speeds  of  15  m  s-1 
results  are  discussed  in  §7. 
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LANGMUIR  CIRCULATIONS  AND  TURBULENT  DISPERSION  OF  OIL  DROPLETS. 


BACKGROUND. 

It  is  generally  agreed  that  the  visible  windrows  that  occur  on  the 
surface  of  oceans  and  lakes  most  frequently  result  from  subsurface  convec¬ 
tive  motions  known  as  Langmuir  circulations.  These  motions,  which  were 
first  described  in  detail  by  l.  Langmuir  (1938),  consist  of  a  series  of 
parallel  counter — rotating  vortices  that  are  aligned  with  the  wind  direc¬ 
tion.  The  vortices  produce  downwelling  and  upwelling  regions  that  lie, 
respectively,  below  lines  of  surface  water  convergence  and  divergence. 
Naturally  occuring  surface  tracer  material,  such  as  seaweed,  foam,  or 
slicks  of  organic  film,  are  swept  into  the  convergence  zones,  creating  a 
visible  surface  pattern  of  bands  or  slicks  parallel  to  the  wind  direc¬ 
tion.  The  associated  surface  streaks  are  readily  detectable  at  the  sur¬ 
face  when  the  wind  speed  exceeds  3m/sec. 

In  this  section,  a  summary  of  field  observations  of  Langmuir  circu¬ 
lations  will  be  given,  followed  by  laboratory  experiments  and  a  descrip¬ 
tion  of  the  physics.  Following  this,  the  physics  of  dispersion  of  parti¬ 
cles  in  a  turbulent  fluid  will  be  described. 

2 . I  Field  Observations 

A  sketch  illustrating  the  general  features  of  Langmuir  circulations 
(which  we  will  abbreviate  by  LC)  is  shown  in  Figure  2.1.  The  surface 
streaks  are  often  very  regularly  spaced;  at  other  times,  they  are  somewhat 
less  regular  but  they  are  always  aligned  close  to  the  wind  direction. 
Since  windrows  lie  above  convergence  zones,  the  spacing  between  two  neigh¬ 
boring  rows  describes  the  horizontal  width  of  a  pair  of  counter-rotating 
cells.  Reported  cell  spacings  (Langmuir  (1938),  Scott  el  al  (1969))  range 
from  2-25  meters  in  lakes,  and  from  2m  up  to  hundreds  of  meters  (at  least 
300m)  in  the  ocean  (Langmuir  (1938), Katz  et  al  (1965),  Ichiye  (1967), 
Gordon  (1970),  Assaf  et  al  (1971)).  The  largest  scales  are  most  easily 
detectable  by  aerial  observation.  Large  scale  LC's  are  associated  with 
smaller  cells  that  are  somewhat  less  well-defined;  thus,  when  large  scales 
appear,  they  often  exist  as  members  of  a  heirachy  of  different  cell  sizes, 
Langmuir  (1938),  Williams  (1965). 
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Figure  2-1.  Sketch  of  Langmuir  circulations  showing  the 
nature  and  scales  of  the  motion. 


Conflicting  information  has  been  given  concerning  factors  that 
determine  cell  width.  A  correlation,  roughly  La4Va  seconds,  between 
cell  width  L  and  wind  speed  Va  has  been  reported  by  some  observers  (Katz 
et  al  (1965),  Maratos  (1971),  Faller  &  Woodcock  (1964),  for  LC  in  the 
ocean.  Others  (Scott  et  al  (1969),  Myer  (1969))  report  work  on  Lake 
George  that  suggests  only  a  slight  correlation  with  wind  speed.  Observa¬ 
tions  on  Lake  George  by  Langmuir  (1938),  Scott  et  al  (1969)  and  Myer 
(1969)),  and  combined  aerial  and  surface  measurements  in  the  ocean  by 
Assaf  et  al  (1971)  suggest  a  correlation  of  row  spacing  with  the  depth  of 
the  seasonal  thermocline. 

It  is  more  difficult  to  determine  the  depth  of  penetration  ;  the 
cells.  The  common  expectation  is  that  the  cell  penetration  is  cor.  ’-able 
to  the  distance  between  windrows.  Some  data  exist  (Langmuir  ;8), 

Scott  et  al  (1969),  Myer  (1969),  Assaf  et  al  (1971))  that  tend  to  .irm 

Langmuir's  conclusion  that  the  cells  can  extend  to  depths  just  above  the 
seasonal  thermocline. 

Observational  information  on  current  speeds  consists  of  disconnected 
parts.  Both  the  cross-wind  (or  sweeping)  component  towards  convergence 
lines,  the  windward  component,  and  the  downwelling  component,  have  been 
measured,  but  not  simultaneously.  From  the  data  reported,  however,  one 
can  estimate  that  the  average  sweeping  speed  is  comparable  to  the  downwel¬ 
ling  speeds.  Data  on  downwelling  speeds  Wj  from  several  sources  have 
been  collected  by  Scott  et  al  (1969),  and  found  to  be  well-correlated  by 
0.0085  Va.  A  very  recent  report  by  Filatov  et  al  (1981)  on  measurements 
made  in  Lake  Ladoga  in  the  Soviet  Union  also  provides  correlation  of  the 
vertical  velocity  near  (the  specific  depth  is  not  reported)  the  water  sur¬ 
face  in  convergence  (downwelling)  zones.  Analysis  of  160  measurements 
with  wind  speeds  in  the  range  3-14  ms~  ,  yields  W<j  =  cjVa  z  where 

Va  is  the  average  wind  speed  10m  above  the  surface,  and  cj  and  C2  appar¬ 

ently  depend  upon  the  density  stratification  of  the  airflow  just  above  the 
surface.  For  neutral  airflow  conditions,  they  cite  the  values 
ci^.bxlO-2 ,  C2a,0.7  (c;  has  dimensions  (cm/s)C2  *),  for  unstable 

airflow  (density  decreasing  upwards)  cia3.7xl0”  ,C2=0.7,  and  for 

2 

stable  airflow  ci=6xl0“  ,  C2!*0.6,  Thus,  all  available  data  indicate 
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that  W<j  increases  with  Va  although  the  rate  is  not  agreed  upon.  Since 
the  average  surface  wind  drift  speed  is  usually  estimated  to  be  around 
0.03Va  to  0.035Va,  and  since  this  is  a  good  measure  of  the  maximum 
mass  transport  speed  in  the  ocean  that  can  be  traced  to  the  wind,  the 
downwelling  speeds  in  LC  are  surprisingly  large.  Upwelling  velocities, 
centered  under  lines  of  surface  divergence,  are  always  reported  to  be 
smaller  than  downwelling  speeds  (implying  an  asymmetry  of  the  cells); 
Langmuir  (1938)  estimated  upwelling  to  be  about  half  as  strong  as  downwel¬ 
ling  motion.  Surface  windward  speeds  are  larger  near  lines  of  surface 
convergence  and  the  effect  is  sufficiently  large  to  be  noted  by  all  obser¬ 
vers  . 

LC  are  readily  observable  for  wind  speeds  exceeding  3m  s-1,  but  this 
is  not  a  critical  condition  and  observations  have  been  made  at  speeds  low¬ 
er  than  0.7  m/s  (Stommel,  1951).  As  Faller  (1981)  has  pointed  out,  the 
commonly  cited  figure  of  3m  s-1  is  open  to  doubt  since  the  observation  of 
streakiness  is  subjective;  the  variable  amount  of  material  available  on 
the  surface  as  tracer,  and  other  factors,  may  play  a  role  in  making  ident¬ 
ification  of  existing  streaks  less  likely  at  lower  winds.  It  is  in  fact 
likely,  as  theory  indicates  (supported  by  laboratory  experiments  de  cribed 
in  the  next  subsection),  that  exceptional  circumsLances  can  occur  in  which 
LC  appear  in  the  absence  of  wind. 

The  time  required  to  form  LC  after  the  wind  begins  to  blow,  or  to 
reorient  the  cells  following  a  shift  in  wind  direction,  is  on  the  order  of 
minutes.  For  example,  Langmuir  (1938)  reported  that  a  large  scale  system 
(spacing  of  hundreds  of  meters)  of  LC  in  the  North  Atlantic  reformed  and 
realigned  with  the  wind  20  minutes  after  a  shift  in  the  wind  direction. 

2 . 2  Laboratory  Experiments 

Laboratory  experiments  on  Langmuir  circulations  have  been  carried 
out  in  wind-wave  tanks.  The  experiments  demonstrated  that  the  phenomenon 

i 

can  be  produced  in  the  laboratory  provided  surface  waves  and  a  sheared 
current,  as  an  applied  wind  stress  necessarily  produces,  are  s imul taneous- 
ly  present.  A  wind  stress  itself  is  not  necessary,  but  the  current  shear 
is  essential;  the  existence  of  a  current  without  waves,  or  waves  without 
current,  is  insufficient  for  the  production  of  LC's. 


8 


The  experiments  of  Faller  (1978)  were  specifically  designed  to  test 
one  of  the  two  mechanisms  (called  CL  I  by  Faller  and  Caponi  (1978))  by 
which  the  Craik  Leibovich  (Craik  and  Leibovich,  1976,  Leibovich  1977a, 
Leibovich  and  Radhakr ishnan,  1977)  theory  produces  Langmuir  circulations. 
The  experiments  of  Faller  and  Caponi  (1978),  like  those  of  Faller  (1978), 
carried  out  in  a  relatively  small  wind-wave  tank  at  the  University  of 
Maryland,  also  had  as  one  primary  purpose  the  testing  of  the  CL  I  mechan¬ 
ism. 

As  will  be  explained  further  in  the  next  subsection,  the  CL  I  mech¬ 
anism  produces  LC  by  direct  wave  forcing  that  can  be  exerted  by  a  coherent 
short -crested  sea.  The  CL  I  theory  shows  that  when  the  wind  stress  is 
applied  in  the  directions  of  wave  propagation  LC  form  with  surface  conver¬ 
gence  lines  at  nodal  lines  of  the  wave  field,  and  divergence  lines  midway 
between;  the  direction  of  circulation  in  the  LC  is  predicted  to  reverse 
when  the  wind  stress  is  applied  in  the  direction  opposite  to  that  of  wave 
propagation.  Regular  short-crested  wave  fields  were  produced  in  these 
experiments  by  specifically  designed  wave-makers,  and  a  wind  stress  was 
also  supplied;  in  all  cases,  the  experiments  confirmed  the  predictions  of 
the  CL  I  mechanism.  The  prediction  of  reversal  of  cell  circulation  was 
confirmed  in  Faller's  (1978)  experiments. 

All  of  these  experiments  were  carried  out  before  the  CL  II  mechan¬ 
ism,  which  produces  LC's  by  an  instability,  (which  can  arise  from  a  com¬ 
pletely  random  wave  field  in  the  presence  of  current  shear  created  by  an 
applied  stress,  to  be  discussed  in  the  next  section)  had  been  fully  de¬ 
scribed.  Nevertheless,  earlier  experiments  by  Faller  (1969)  showed  that 
LC  could  be  produced  in  the  laboratory  by  random  wind-waves.  Furthermore, 
the  scalings  of  both  CL  I  and  CL  II  theories  are  the  same,  and,  as  it 
turns  out  (Leibovich  and  Paolucci  (1980)),  the  structure  and  strengths  of 
the  motions  in  LC  cells  produced  by  the  two  mechanisms  are  very  similar. 

2.3.  Physics  of  Langmuir  Circulations 

Although  a  variety  of  ideas  has  been  suggested  to  explain  Langmuir 
circulations  (see  Scott  et  al,  1969,  for  example,  for  a  discussion  of  most 
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of  these  early  specifications),  it  j  now  generally  agreed  that  all  of 
these  processes,  except  those  involving  interactions  between  waves  and 
currents,  are  either  irrelevant  (cannot  account  for  LC)  or  not  necessary 
but  may  modify  the  phenomenon  (i.e.  thermal  convection  can  strenghten  LC, 
but  since  LC  are  observed  in  thermally  stable  conditions,  thermal  insta¬ 
bility  is  not  an  underlying  causal  mechanism  for  LC).  Furthermore,  it  now 
seems  clear  (Leibovich,  1980)  that  all  aspects  of  wave  and  current  inter¬ 
actions  relevant  to  LC  (for  example,  Garrett's  (1976)  theory),  may  be 
incorporated  in  the  CL  (Craik-Leibovich)  theory,  which  will  now  be 
described . 

If  attention  is  focussed  on  phenomena  occuring  on  a  time  scale 
larger  than  a  typical  period  of  a  surface  gravity  wave,  it  is  appropriate 
to  filter  out  the  wave-generated  "hash"  by  considering  only  quantities 
averaged  over  several  wave  periods.  The  averaged  momentum  equations  con¬ 
tain  the  rectified  effects  of  the  surface  waves  only  through  a  modifica¬ 
tion  of  the  pressure,  and  through  an  apparent  force  term  which  represents 
a  direct  interaction  between  the  waves  and  current.  The  first  effect  is 
of  not  particular  consequence  since  it  only  alters  the  mean  pressure  by  a 
small  amount.  The  second  effect,  however,  is  crucial.  Wave-current 
interactions  arise  as  a  vortex  force,  (Leibovich,  1977b)  which  may  be 
represented  as 


=  u  x  curl  u8 


(2-1) 


where  u8  is  the  particle  drift  velocity  vector 
caused  by  progressive  water  waves  and  u  =  (u,v,w) 
city  vector  measured  at  the  point  x,  y,  z  at  time 
of  the  waves  is  carried  by  the  Stokes  drift.) 


(Stokes  drift,  see  §4) 
is  the  mean  water  velo- 
t.  (Thus,  the  momentum 
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Suppose  we  consider,  for  simplicity  of  explanation,  a  case  where  u 
and  us  are  both  in  the  same  direction  (say  the  x-direction)  and  both 
decrease  with  depth  (where  z  is  measured  vertically  upwards  from  the  mean 
free  surface)  but  do  not  vary  horizontally.  This  situation  is  the  sim¬ 
plest  intuitive  model  of  a  wind-generated  sea;  the  average  motion  is  in 
the  wind  direction  and  varies  only  with  depth,  while  the  wave  field 
(statistically)  does  not  vary  in  the  horizontal  (thus  leading  to  a  Stokes 
drift  varying  only  with  depth.  Then  Fv  has  a  component  only  in  the 
z-direction,  (upwards,  acting  like  a  buoyant  force),  and  therefore  has  an 
effect  analogous  to  density  stratification.  Provided  the  symmetry  of  a 
perfectly  layered  configuration  is  attained,  the  motion  can  continue  to 
maintain  its  layered  structure.  If,  however,  there  is  a  slight  deviation 
e.g.,  if  any  horizontal  variations  in  either  u  or  u8  develop,  the 
vortex  force  causes  accelerations  perpendicular  to  the  current  direction 
and  Langmuir  circulations  must  develop.  This  we  show  below,  assuming  for 
simplicity  that  there  are  no  variations  in  the  direction  of  the  wind. 

This  is  most  easily  seen  when  the  Stokes  drift  varies  horizontally; 
in  that  event  LC  are  created  by  direct  wave  forcing:  this  is  the  CL  I 
mechanism  referred  to  earlier.  For  short  time  intervals  at  least,  the 
structure  of  wind  waves  is  typically  short -crested,  with  the  waves  appear¬ 
ing  to  be  due  to  two  coherent  wavetrains  propagating  at  equal  but  opposite 
angles  to  the  wind  direction.  This,  if  maintained  for  a  time  long  enough 
to  obtain  meaningful  averages,  leads  to  a  spanwise  periodic  Stokes  drift, 
(see  Figure  2.2),  so  that  9us/3y  *  0.  This  leads  to  an  acceleration  in 
the  spanwise  (y)  direction  perpendicular  to  the  current  of  amount 

-(3ug/3y) (3u/3z) 


Since  3ug/3y  alternates  in  sign  the  water  moves,  from  both  sides, 


towards  those  planes  where  us  is  a  minimum.  There,  by  continuity,  water 
must  sink. 

The  CL  I  mechanism  requires  that  the  wave  structure  remain  coherent 
for  time  scales  comparable  with  the  lifetime  of  a  Langmuir  cell.  While 
this  may  in  fact  be  the  case  under  certain  circumstances,  it  is  not  likely 
to  be  common.  The  second,  or  CL  II,  (Craik,  1977,  Leibovich  1977a  and 
Paolucci  1980a,  1980b,  1981),  mechanism  is  an  instability  and,  by  con¬ 
trast,  is  always  a  possibility.  Suppose  u8  is  horizontally  uniform 
(consistent  with  statistically  homogeneous  waves),  but  that  a  small  span- 
wise  irregularity  is  present  in  an  otherwise  horizontally  uniform  cur¬ 
rent.  Such  irregularities,  of  course,  are  unavoidable  and  everywhere 
present  in  nature.  In  such  a  case  3u/3y  *  0,  and  the  vortex  force 

Fv  causes  an  acceleration  towards  the  planes  of  maximum  u,  where  by 
continuity  the  fluid  must  sink  (Figure  2-3). 

It  is  essential  to  know  whether  the  initial  irregularity  grows.  If 
not,  then  there  is  no  reason  for  the  horizontal  motions  to  grow,  and  all 
that  we  conclude  is  that  these  small  irregularities  in  the  wind  directed 
velocity  component  are  accompanied  by  comparably  small  irregularities  in 
the  other  two  velocity  components.  If  a  feedback  mechanism  exists  permit¬ 
ting  this  system  to  grow,  then  coherent  and  measureable  convection  pat¬ 
terns  can  develop.  Conservation  of  x-momentum  for  a  vanishingly  thin  slab 
of  fluid  surrounding  the  convergence  plane  shows  that  as  the  fluid  sinks 
the  x-velocity  must  increase,  assuming,  that  3u/3z  >  0  (as  one  expects  in 
an  oceanic  current).  Thus,  provided  this  increase  of  x-speed  is  not  elim¬ 
inated  by  frictional  effects,  the  slight  acceleration  of  fluid  leading  to 
the  initial  current  irregularity  will  force  side  motions  due  to  vortex 
force  and  this  will  in  turn  amplify  the  irregularity.  We  therefore  have 
an  instability  of  the  system  which  has  the  potential  to  develop  into  LC. 
The  analysis  in  the  literature  (in  particular,  Leibovich  and  Paolucci 
1981),  shows  that  in  typical  wind-generated  seas,  frictional  forces  almost 
never  are  strong  enough  to  overcome  this  feedback  process. 

In  chapter  4,  we  will  present  the  theory  in  mathematical  form.  The 
CL  I  and  CL  II  mechanisms  are  based  upon  the  same  scaling  assumption  and 
therefore  the  general  features  of  the  convective  cells  resulting  from  both 
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Figure  2-3.  Sketch  illustrating  the  generation  of  Langmuir 
circulations  by  the  instability  mechanism.  The 
Stokes  drift  need  not  vary  in  the  spanwise  direc¬ 
tion. 
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the  direct  and  the  instability  mechanisms  are  similiat  (Leibovich  and 
Radhakr ishnan,  1977;  Craik,  1977;  Leibovich  and  Paolucci,  1980a).  Since 
the  instability  mechanism  does  not  depend  upon  any  special  surface  wave 
organization,  it  would  seem  to  be  the  more  attractive  possibility  for 
Langmuir  circulations  and  this  is  the  mechanism  utilized  in  the  present 
work.  It  should  be  noted,  however,  that  the  direct  forcing  mechanism 
always  acts,  but  may  be  able  to  exert  a  coherent  effect  only  for  short 
periods  of  time.  While  this  short-time  duration  limits  its  ability  to 
produce  a  long-lived  Langmuir  cell  structure,  the  direct  forcing  can  be 
expected  to  lead  to  initial  perturbations  that  can  be  subsequently 
amplified  into  long-lived  Langmuir  circulations  by  the  instability 
mechanism. 

2.4  Description  of  the  Physics  of  Turbulent  Dispersion  of  Oil  Droplets. 

When  there  is  relative  motion  between  an  oil  droplet  and  the  sur¬ 
rounding  sea  water,  a  drag  is  produced.  This  drag  may  be  steady,  un¬ 
steady,  or  quasi-steady  (which  we  will  explain  in  a  moment),  depending  on 
the  circumstances.  The  drag,  of  course,  results  from  the  fact  that  the 
sea  water  cannot  slip  where  it  meets  the  surface  of  the  oil  droplet;  if 
there  are  impurities  on  the  interface  between  the  seawater  and  the  oil, 
the  droplet  skin  may  act  like  a  rigid  surface  to  a  greater  or  lesser  ex¬ 
tent.  Otherwise,  there  may  be  a  circulation  within  the  oil  droplet  so 
that  there  is  art  effective  slip  at  the  surface,  although  the  velocity  of 
circulation  within  the  droplet  will  be  relatively  small  in  any  event, 
since  the  oil  is  so  much  more  viscous  than  the  sea  water.  The  droplet  may 
deform,  depending  on  the  size,  velocity  and  surface  tension,  combined  in  a 
dimensionless  group  called  the  Weber  number.  Under  more  extreme  condi¬ 
tions  it  may  be  torn  apart.  It  may  interact  with  other  droplets,  depend¬ 
ing  on  the  concentration.  Depending  on  the  Reynolds  number  of  the  rela¬ 
tive  flow  around  the  droplet,  the  drag  may  be  a  non-linear  function  of  the 
relative  velocity,  and  other  non-linear  effects  may  become  important,  such 
as  lateral  forces  resulting  from  spin  (the  Magnus  effect),  or  even  from 
the  unsteadiness  of  the  wake,  which  may  result  in  the  drag  not  being  in 
the  instantaneous  direction  of  the  relative  motion.  If  the  oil  has  lost 
its  fluidity,  and  forms  into  non-spherical  oil  particles  rather  than 
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droplets,  slow  rotation  of  the  particles  will  result  in  unsteady  lateral 
forces,  as  the  particle  deflects  the  oncoming  flow  now  to  one  side,  now  to 
the  other.  At  higher  values  of  the  relative  Reynolds  number  unsteady  sep¬ 
aration  may  occur,  leading  to  unsteady  lateral  forces.  A  material  point 
(of  water)  in  a  turbulent  flow  will  slowly  wander  from  its  original  neigh¬ 
borhood;  that  is,  in  a  turbulent  flow  under  waves,  the  original  neighbor¬ 
hood  will  follow  an  orbit  described  by  the  potential  motion  driven  by  the 
surface  waves.  A  material  point  will  begin  by  tracing  one  of  these  or¬ 
bits,  but  will  gradually  wander  away  from  it.  An  oil  droplet  will  do  a 
similar  thing,  but  will  also  wander  away  from  the  material  point,  since 
there  is  a  relative  motion  between  it  and  the  sea  water.  Hence,  the  dif¬ 
fusion  of  the  oil  droplet  will  be  different  from  the  diffusion  of  the 
material  point.  In  addition,  the  mean  density  of  the  oil  droplet  is  dif¬ 
ferent  from  that  of  the  sea  water.  As  a  result,  the  droplet  will  rise  (in 
general)  through  the  water.  In  rising,  it  leaves  the  parcel  of  sea  water 
with  which  it  started,  and  hence  the  velocity  history  following  the 
wandering  droplet  will  in  general  lose  correlation  faster  than  that  of  a 
wandering  material  point.  Finally,  if  we  include  a  mean  flow  (such  as  the 
Langmuir  circulation),  we  envision  a  droplet  with  a  basic  motion  consist¬ 
ing  of  a  superposition  of  a  very  small  steady  rise,  a  relatively  rapid 
orbital  motion,  the  Langmuir  circulation  (stronger  than  the  steady  rise), 
and  superposed  on  this  an  unsteady  wandering  resulting  from  both  the  un¬ 
steady  turbulent  velocity  and  the  other  possible  sources  of  unsteadiness 
inherent  in  the  droplet's  own  wake,  and  in  the  interactions  with  other 
droplets . 

This  is  clearly  an  extremely  complicated  situation.  No  approach 
exists  which  is  capable  of  handling  all  these  phenomena  simultaneously. 
Fortunately,  it  is  possible  to  sort  realistic  scenarios  into  groups  in 
which  one  or  more  of  the  complicating  factors  are  not  present,  and  obtain 
a  solution.  Scenarios  which  cannot  be  handled  in  this  way  represent  rela¬ 
tively  rare  occurrences. 

The  first  difficulty  to  be  eliminated  is  interaction.  There  is  no 
method  for  dealing  with  interacting  droplets.  For  spherical  or  non-spher- 
ical  droplets,  interaction  is  negligible  for  volume  fractions  (of  circum- 
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scribed  spheres)  below  0.003  (Luraley,  1978).  To  get  a  physical  feel  for 
this  number,  for  oil  droplets  one  millimeter  in  diameter,  this  volume 
fraction  corresponds  to  six  droplets  per  cubic  centimeter.  While  it  is 
clear  that  many  situations  of  interest  may  have  volume  fractions  higher 
than  this,  it  is  also  clear  that  many  situations  of  interest,  in  particu¬ 
lar  the  transport  and  dispersion  of  droplets  by  turbulence  in  Langmuir 
cells,  are  likely  to  satisfy  this  limitation.  We  will  probably  violate 
this  limitation  only  near  the  windrow  area,  where  we  expect  the  droplets 
to  coalesce  in  any  event,  so  that  a  droplet  description  is  no  longer 
valid.  In  everything  that  follows,  we  will  assume  negligible  interac¬ 
tion.  This  does  not  mean  that  there  is  no  interaction;  only  that  the 
interactions  are  sufficiently  infrequent  for  the  presence  of  the  droplets 
to  make  a  less  than  one  percent  contribution  to  the  viscosity  of  the 
suspension.  Hence,  it  would  be  legitimate  to  use  a  non-interactingmodel 
to  compute  the  rate  of  coalescence  of  droplets  at  low  concentrations,  for 
example . 

The  next  difficulty  is  tfyat  of  the  unsteady  wakes.  There  is 
no  known  method  for  treating  vortex  shedding,  unsteady  separation,  Magnus 
effects,  and  the  like.  Fortunately,  these  are  all  Reynolds-number  ef¬ 
fects,  and  extensive  information  exists  on  the  relative  Reynolds  number  at 
which  the  wake  becomes  unsteady  (see  Van  Dyke,  1975).  This  appears  to  be 
approximately  60.  Hence,  as  long  as  we  consider  only  particles  for  which 
the  Relative  Reynolds  number  is  less  than  60,  we  can  eliminate  most  of  the 
unsteady  effects.  The  only  exception  is  the  fluctuating  sideways  force 
resulting  from  the  Magnus  effect,  which  we  wi.ll  discuss  in  a  moment.  In 
section  5.  we  will  discuss  in  detail  the  implications  of  this  limiting 
relative  Reynolds  number  for  droplet  size;  the  relative  Reynolds  number  is 
due,  of  course,  to  both  droplet  inertia  and  to  the  terminal  velocity  re¬ 
sulting  from  buoyancy.  Different  density/diameter  combinations  result  in 
different  situations.  In  section  5  we  will  conclude  that  droplets  of 
roughly  1.5  mm  in  diameter  can  be  considered  without  exceeding  the  Rey¬ 
nolds  number  limitation.  Droplets  of  this  size  certainly  include  many  of 
those  of  interest. 

So  far  as  the  Magnus  effect  is  concerned,  a  droplet  in  a  shear  will 
try  to  spin  with  the  local  vorticity.  To  the  extent  that  the  flow  field 
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around  the  droplet  is  non-linear  (and  it  certainly  will  be  if  the  relat ive 
Reynolds  number  is  not  small  compared  to  unity),  if  there  is  a  relative 
velocity  there  will  be  developed  a  force  normal  to  the  direction  of  the 
relative  velocity.  We  may  estimate  the  order  of  magnitude  of  this  force 
as 

Pain'd  (2-3) 

where  pw  is  the  density  of  the  water,  AU  is  the  relative  velocity,  and  T 
is  the  circulation  about  the  droplet.  This  is  the  standard  expression 
used  for  the  lift  on  wings  and  other  bodies  which  induce  a  circulation. 
We  may  estimate  T  as  nd  u/X,  where  u  is  the  rms  fluctuating  velocity  of 
the  turbulence,  and  X  is  the  Taylor  microscale.  This  assumes  that  the 

droplet  will  rotate  with  the  local  shear.  Equivalently,  we  may  use  the 

expression  3 . 8  7  A  /  Rjj.  1/2  for  \  where  l  is  the  length  of  the  energy-con¬ 
taining  eddies,  and  Re  is  the  Reynolds  number  based  upon  u  and  l.  If  we 

take  the  ratio  of  the  lift  force  (  say,  Fp)  to  the  drag  force  (say, 
FD)  we  obtain 


Fl/Fd  *  3(pw/pQ)(ua/X)/[l  +  (3/16)R]  (2-4) 

where  R  is  the  relative  Reynolds  number,  a  is  the  droplet  time  scale,  and 
Pq  is  the  density  of  the  oil  droplet.  The  substantive  part  of 
expression  (2-4)  is  the  ratio  of  the  droplet  time  scale  to  the  time  scale 
of  the  dissipative  eddies.  This  is  the  same  time  scale  ratio  that  will  be 
estimated  later  in  another  context,  and  it  is  a  small  number.  We  will 

Q  , 

find  that  a  is  no  larger  than  3x10"  sec,  that  a  reasonable  estimate  for  u 
is  1x10"  m/s,  and  that  Z  is  roughly  1  m,  while  Rjj  is  about  1x10“  , 
giving  X/u  *  3.87  sec.  Hence,  we  may  conclude  that  the  fluctuating  lift 
forces  g_.ierated  by  the  spin  are  a  very  small  fraction  of  the  drag  forces, 
and  hence  that  the  dispersion  produced  by  these  forces  is  a  similar 
fraction  of  the  dispersion  produced  by  the  drag  forces. 

One  final  source  of  unsteadiness  and  lateral  forces  cannot  be  com¬ 
pletely  eliminated.  If  the  oil  particle  is  non-spherical  (we  cannot  now 
say  droplet),  a  situation  that  may  arise  depending  on  the  age  and  condi¬ 
tion  of  the  oil,  the  drag  will  not  be  steady  and  will  not  be  instantan¬ 
eously  in  the  direction  of  the  relative  motion.  Even  in  the  absence  of 
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turbulence,  the  oil  particle  will  slowly  rotate  as  it  rises,  and  as  it 
rotates  the  direction  and  magnitude  of  the  drag  will  fluctuate  as  the  oil 
particle  presents  different  aspects  to  the  oncoming  flow.  The  lateral 
force  resulting  from  this  effect  is  of  the  same  form  as  that  due  to  spin, 
discussed  above,  but  here  the  circulation  is  caused  by  the  lateral  de¬ 
flection  of  the  flow  by  the  orientation  of  the  oil  particle.  We  may  guess 
that  the  maximum  velocity  difference  across  the  oil  particle  that  can  be 
induced  in  this  manner  will  be  no  larger  than  the  oncoming  relative  velo¬ 
city.  With  this  assumption,  the  lateral-to-drag  force  ratio  has  the  same 

form  as  in  equation  (2-2),  but  the  important  part  is  now  (AUa/d).  The 

-2 

maximum  relative  velocity  will  be  estimated  as  no  more  than  roughly  1x10 

3  2 

m/s,  with  d  =  1x10"  m,  so  that  the  ratio  is  about  3x10"  ;  hence,  the 

lateral  forces  resulting  from  lateral  turbulent  gusts  are  much  larger  than 

the  lateral  forces  resulting  from  asymmetrical  orientation  of  non-spheri- 

cal  oil  particles.  In  what  folLows,  we  will  refer  generally  to  droplets, 

implying  sphericity,  but  in  fact  what  we  say  will  be  applicable  to  non- 

spherical  oil  particles.  The  appropriate  terminal  velocity,  of  course, 

cannot  be  determined  from  the  expressions  for  spherical  droplets,  but  must 

in  general  be  determined  experimentally. 

The  next  problem  that  we  must  eliminate  has  to  do  with  unsteady  ef¬ 
fects  in  the  drag  (other  than  those  discussed  above).  That  is,  our  prob¬ 
lem  will  be  made  much  easier  if  the  drag  may  be  assumed  to  be  a  quasi¬ 
steady  function  of  the  relative  velocity.  Unsteady  effects  of  this  type 
result  from  the  time  necessary  to  establish  a  viscous  region  around  the 
particle  in  equilibrium  with  the  oncoming  flow.  The  thicker  the  viscous 
region,  generally  speaking,  the  longer  it  takes  to  establish  it  when  there 
is  a  change  in  the  oncoming  relative  velocity.  Hence,  for  droplets  in  the 
Stokes  regime  (relative  Reynolds  numbers  small  compared  to  unity)  this  can 
be  a  problem,  and  for  density  ratios  close  to  unity  the  additional  forces 
resulting  from  this  effect  can  be  non-negl igible .  This  type  of  force  is 
referred  to  as  a  Basset  force,  from  the  investigator  who  first  discussed 
the  form  of  this  term  in  the  equation.  Fortunately,  the  vast  majority  of 
the  oil  droplets  in  which  we  are  interested,  the  larger  ones  which  carry 
most  of  the  oil  volume,  have  a  relative  Reynolds  number  well  above  the 
Stokes  range.  For  these  droplets,  the  viscous  region  is  relatively  thin, 
and  consequently  the  time  required  to  establish  it  after  a  change  in  the 
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relative  velocity  is  negligibly  short.  In  Lumley  (1978)  estimates  are 
made  of  the  relative  time  scales;  we  do  not  reproduce  them  here,  but 
simply  state  that  for  al’  the  droplets  of  interest  to  us,  we  may  consider 
that  the  drag  is  a  quasi-steady  function  of  the  relative  velocity,  so  that 
the  drag  force  may  always  be  taken  to  be  in  equilibrium  with  the  relative 
velocity. 

The  next  complication  is  that  of  non-linear  drag.  If  the  relative 
Reynolds  number  is  as  large  as  sixty,  it  is  clear  that  the  drag  force  is 
not  linear  in  the  relative  velocity.  This  could  be  a  very  complicating 
factor;  if  the  droplet  rise  were  opposed  by  a  turbulent  gust,  the  downward 
force  would  be  greater  than  the  corresponding  upward  force  if  a  similar 
gust  were  to  aid  the  rise.  Hence,  the  terminal  velocity  would  be  reduced 
in  a  turbulent  field,  and  would  be  a  function  of  the  turbulence  paramet¬ 
ers,  as  well  as  the  droplet  inertia.  In  general,  the  droplet  dispersion 
would  be  dominated  by  the  more  intense  turbulent  velocity  fluctuations, 
and  prediction  of  the  dispersion  would  become  very  difficult.  However,  if 
we  could  linearize  around  the  rise  at  the  terminal  velocity,  all  these 
problems  would  be  eliminated.  This  would  automatically  eliminate  the 
change  in  the  terminal  velocity  due  to  the  turbulence,  and  would  produce  a 
fluctuating  drag  force  which  was  a  linear  function  of  the  turbulent  fluc¬ 
tuations,  making  the  prediction  of  the  dispersion  relatively  simple. 
Fortunately,  the  conditions  for  this  linearization  are  simple  to  obtain, 
and  they  result  in  a  fairly  loose  restriction  on  the  fluctuating  relative 
Reynolds  number.  We  find  that  as  long  as  the  fluctuating  relative  Rey¬ 
nolds  number  is  less  than  12.5,  the  first  nonlinear  term  neglected  is  no 
more  than  20%  of  the  linear  terra;  hence,  in  the  worst  case  we  will  make  no 
more  than  a  20%  error  in  the  dispersion.  Considering  the  other  uncertain¬ 
ties  in  this  analysis,  a  20%  error  is  certainly  acceptable.  This  error, 
in  fact,  will  in  general  be  much  smaller. 

The  prediction  of  the  fluctuating  relative  Reynolds  number  is  diffi¬ 
cult,  because  the  seawater  velocity  at  the  droplet  location  depends  on  the 
droplet  motion.  This  is  a  classical  problem  which  has  not  been  completely 
resolved.  However,  two  limiting  cases  are  relatively  simple  to  under¬ 
stand:  the  case  in  which  the  droplet  essentially  follows  the  motion  of  the 
sea  water  perfectly,  called  the  Lagrangian  case,  and  the  case  in  which  the 
droplet  nearly  ignores  the  motion  of  the  sea  water,  rising  through  it  rap- 
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idly  (relatively).  Of  coarse,  both  cases  are  approximations,  and  the 
truth  lies  somewhere  between;  however,  consideration  of  the  two  cases  al¬ 
lows  bracketing  estimates  to  be  made  of  the  fluctuating  relative  Reynolds 
number.  Proceeding  in  this  way,  we  will  find  that  the  majority  of  the 
droplets  in  which  we  are  interested  are  much  closer  to  the  Eulerian  limit 
than  they  are  to  the  Lagrangian  limit.  We  easily  obtain  in  this  way  a 
limitation  on  droplet  diameter  to  assure  the  appropriate  limitation  on  the 
fluctuating  relative  Reynolds  number,  and  find  that  we  can  consider  drop¬ 
lets  up  to  approximate  1.5  millimeters  in  diameter. 

Finally,  we  must  consider  the  fact  that  the  wandering  droplet  rises 
out  of  the  turbulent  eddy  with  which  it  begins,  and  the  fluctuating  drop¬ 
let  velocity  loses  correlation  faster  than  the  fluctuating  velocity  of  a 
material  point  which  began  with  it.  This  is  called  the  crossing-traject¬ 
ories  effect,  and  has  been  much  studied  by  Csanady  and  Yudine  (for  refer¬ 
ences,  see  Lumley,  1978).  Essentially,  the  matter  is  resolved  in  an  ap¬ 
proximate  manner  by  noting  that  in  the  two  limiting  cases  discussed  above 
the  dispersion  takes  on  a  relatively  simple  form.  Essentially,  an  inter¬ 
polation  form  between  these  two  limiting  cases  is  constructed,  assuming 
only  that  the  velocity  correlations  at  all  intermediate  values  of  the 
parameters  are  similar  in  form,  and  that  iso-correlation  contours  are  el¬ 
liptical.  The  resulting  forms  are  particularly  simple,  and  depend  only  on 
the  ratio  of  the  terminal  velocity  to  the  turbulent  fluctuating  velocity. 
The  net  result  is  a  reduction  of  the  dispersion  of  the  droplet  relative  to 
the  dispersion  of  a  material  point,  and  a  difference  of  the  dispersion  in 
the  vertical  direction  relative  to  that  in  the  horizontal  direction. 

Finally,  we  must  consider  the  dispersion  as  a  function  of  time.  As 
a  droplet  begins  to  be  dispersed  by  the  turbulent  fluctuations,  the  dis¬ 
persion  changes  as  a  function  of  time,  growing  at  first  parabolically,  and 
later  linearly  in  time.  The  exact  form  of  the  growth  depends  on  the  form 
of  the  correlation  functions  that  is  assumed;  it  is  possible  to  make 
simple  assumptions  that  are  more  than  adequate  for  this  prediction.  How¬ 
ever,  if  this  transient  period  is  a  relatively  small  part  of  the  total 
time  over  which  the  droplets  are  dispersed  by  the  turbulence,  the  tran¬ 
sient  period  may  be  ignored,  and  the  asymptotic  value  of  the  diffusivity 
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used;  the  error  introduced  in  this  way  is  negligible.  In  our  case,  for¬ 
tunately,  this  is  the  case,  since  the  circulation  in  the  Langmuir  cells  is 
relatively  slow,  and  the  transient  period  of  the  dispersion  corresponds  to 
a  small  fraction  of  the  cell  turnover  time. 

In  this  way  we  have  finally  reduced  the  droplet  dispersion  problem 
to  the  dispersion  of  spherical,  non-interacting  droplets  with  constant 
diffusivities,  different  in  different  directions,  and  under  different  con¬ 
ditions.  Such  a  situation  may  be  described  by  a  diffusion  equation  in 
which  the  diffusivities  are  functions  of  position  and  direction. 
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3. 


PLAN  OF  WORK 


The  primary  purpose  of  this  work  is  the  determination  of  subsurface 
oil  dispersion  created  by  Langmuir  circulations  and  smaller  scale  turbu¬ 
lence.  In  this  section,  the  major  assumptions  adopted  in  the  work  are 
recorded,  the  sequence  of  steps  followed  is  outlined,  and  specific  choices 
are  made  for  the  Langmuir  circulations  to  be  computed. 

Dispersion  results  from  the  random  displacement  of  material  caused 
by  the  fluid  turbulence.  A  description  requires  information  about  the 
turbulence  as  well  as  mean  fluid  motions  and  both  turbulence  and  the  large 
scale  convective  motions  are  extremely  complicated.  According  to  the 
theory  of  Langmuir  circulations  that  we  employ  here,  convective  activity 
arises  as  an  instability  whose  detailed  evolution  depends  upon  the  wind 
stress,  fetch,  ambient  currents,  sea  state,  water  density  stratification, 
water  depth,  turbulent  momentum  transport,  initial  perturbations,  time  of 
observation,  and  the  rotation  of  the  earth.  Some  of  these  effects  have 
been  studied,  but  many  have  not.  Even  if  all  contributing  factors  were 
well  understood,  a  sample  of  Langmuir  circulation  observations  taken  under 
identical  conditions  of  the  principal  environmental  observables  (wind 
speed  and  sea  state)  is  expected  to  yield  random  results  and  only  statis¬ 
tical  measures  of  the  observations  will  be  reproducible.  This  is  because 
the  theory  is  an  instability  theory,  because  the  flow  details  depend  upon 
prior  history,  and  because  many  factors  other  than  the  principal  obser¬ 
vables  influence  the  motion.  This  situation  is  characteristic  of  all  geo¬ 
physical  fluid  dynamics  problems.  One  has  a  reasonable  right  to  expect 
from  a  theory  the  significant  trends;  it  is  not  reasonable  to  expect  an 
accurate  prediction  of  a  real  geophysical  event  when,  as  is  invariably  the 
case,  only  part  of  the  data  required  to  specify  the  event  is  known. 

In  the  present  work,  we  assume  that  the  ocean  is  infinitely  deep, 
and  we  ignore  density  stratification  and  the  earth's  rotation.  The  motion 
is  assumed  to  start  from  rest  under  the  action  of  an  applied  stress,  and 
therefore  it  is  assumed  that  there  are  no  ambient  currents,  or  if  such 
currents  are  present,  the  associated  vorticity  may  be  ignored.  Although 
we  use  a  semi-empirical  model  to  relate  the  sea  state  to  the  wind  speed 
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and  fetch,  we  do  not  take  into  account  the  time  required  for  the  surface 
wave  field  to  develop. 

We  believe  that  Coriolis  accelerations  have  a  major  effect  on  Lang¬ 
muir  circulations  when  the  wind  blows  for  a  time  exceeding  a  quarter  of  a 
pendulum  day.  Although  density  stratification  strongly  modifies  the  rate 
of  penetration  of  Langmuir  circulations  in  deep  water,  it  is  the  Coriolis 
acceleration  which  probably  determines  the  maximum  depth  of  penetration. 
Unfortunately,  this  effect  is  not  yet  understood.  Without  a  Coriolis 
force,  the  motion  in  infinitely  deep  water  has  no  steady  limit,  and  some 
criterion  must  be  used  to  describe  when  to  terminate  a  calculation.  We 
use  arbitrary  specified  (maximum)  penetration  depths  to  help  organize  our 
calculations. 

The  first  step  in  our  analysis  is  the  calculation  of  several  Lang¬ 
muir  circulation  flow  patterns  thought  to  be  representative  of  a  range  of 
environmental  conditions  of  practical  interest;  the  specific  parameter 
choices  are  described  in  §3.1.  Turbulence  is  represented  in  the  Langmuir 
circulations  by  a  constant  eddy  viscosity,  and  the  choices  made  are  also 
described  in  §3.1.  Calculation  of  all  three  velocity  components  requires 
the  solution  of  a  fully  nonlinear  problem.  This  is  carried  out  by  a 
numerical  scheme  to  be  described  in  §4,  and  the  results  are  stored  on  mag¬ 
netic  tape  for  further  use. 

It  is  assumed  that  the  presence  of  oil  has  no  dynamical  conse¬ 
quences,  and  that  the  Langmuir  circulation  fields  are  therefore  not 
affected  by  the  presence  of  oil.  While  this  assumption  is  hard  to  dispute 
as  far  as  subsurface  oil,  which  typically  has  very  low  volume  concentra¬ 
tions,  is  concerned,  the  presence  of  surface  oil  clearly  has  an  effect  on 
surface  waves  and  perhaps  therefore  on  LC.  Nevertheless,  the  assumption 
is  expected  to  be  sound  under  most  circumstances,  since  thin  surface  oil 
films  primarily  damp  short  (capillary)  surface  waves.  Waves  in  the 
gravity  wave  regime  are  responsible  for  LC,  and  these  are  not  expected  to 
be  substantially  affected  by  surface  oil.  We  have  not,  however,  analyzed 
this  question  in  detail,  and  it  is  a  point  that  requires  and  deserves  fur¬ 
ther  study. 
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Each  of  the  computed  Langmuir  circulation  fields  is  regarded  as  a 
possible  outcome  of  an  actual  observation.  Furthermore,  in  order  to  make 
the  dispersion  calculation  tractable,  we  assume  that  these  Langmuir  circu¬ 
lation  fields  are  possible  steady  motions,  even  though  they  are  stages  in 
the  time  development  of  an  inherently  unsteady  mathematical  problem. 
Since  the  motions  generated  are  believed  to  scale  properly  with  the  pri¬ 
mary  environmental  parameters,  we  believe  that  dispersion  calculations 
based  upon  these  assumptions  will  be  at  least  qualitatively  meaningful. 
As  we  have  already  stated,  this  is  all  one  can  usually  expect  in  geophys¬ 
ical  problems. 

After  generating  the  water  motions,  we  consider  the  trajectory  of  a 
single  buoyant  particle  released  at  an  arbitrary  point  in  any  of  our  Lang¬ 
muir  circulation  fields.  The  key  assumption  made  is  that  the  velocity 
field  into  which  the  particle  is  placed  is  completely  time  independent 
(e.g.,  there  are  no  turbulent  fluctuations).  This  idealized  problem, 
which  follows  a  similiar  study  by  Stommel  (1949),  reveals  the  presence  of 
subsurface  zones  in  which  buoyant  particles  can  be  trapped.  If  a  particle 
is  placed  within  a  closed  boundary  defining  one  of  these  zones,  it  will 
move  in  a  closed  orbit  lying  inside  the  boundary  of  the  zone.  The  size 
(depth  and  volume)  of  these  zones  depends  upon  the  buoyancy  of  the  parti¬ 
cle  and  upon  the  intensity  and  length  scale  of  the  LC  motion,  but  for 
typical  values  expected  in  oil-spill  situations,  these  regions  are  exten¬ 
sive.  For  each  of  the  LC  fields  generated,  these  zones  are  mapped  for 
three  values  of  particle  buoyant  terminal  velocity,  and  stored  on  computer 
tape. 

The  generation  of  LC  fields  and  subsurface  trapping  zones  (which  we 
refer  to  as  SRZ,  or  "Stommel  retention  zones")  is  described  in  §3.1  and  in 
§4. 

How  does  a  particle  get  into  a  SRZ?  In  the  idealized  problem,  a 
particle  is  either  in  or  out;  it  cannot  cross  the  closed  SRZ  boundary.  Of 
course,  if  a  particle  can  get  in,  it  must  also  be  possible  to  get  out,  and 
so  the  trapping  is  not  complete.  The  idealized  analysis  omits  velocity 
fluctuations,  and  crossing  of  the  boundary  is  expected  to  be  possible  as  a 
result  of  the  turbulent  fluctuations,  inevitably  present. 


25 


The  investigation  of  the  appropriate  model  for  transport  due  to  tur¬ 
bulence  is  the  subject  of  §5.  This  model  must  allow  the  calculation  of 
the  distribution  of  oil  particles  resulting  from  the  combined  effects  of 
Langmuir  cell  mean  motions  and  superposed  random  velocity  fluctuations. 
The  model  adopted  as  a  consequence  of  the  considerations  of  §5  is  based 
upon  eddy  diffusivity,  and  describes  the  volume  fraction  of  oil  in  terms 
of  convection  by  the  mean  (LC)  velocity  field,  buoyant  rise,  and  diffusion 
due  to  turbulence.  The  principal  matter  to  be  determined  in  such  a  model 
is  the  relationship  between  eddy  diffusivity  and  the  turbulence  character¬ 
istics.  In  the  course  of  §5,  influence  of  wave  accelerations,  and  the 
influence  of  oil  particle  inertia  terras  upon  the  rate  of  approach  of  the 
turbulent  diffusivity  to  its  asymptotic  values  are  evaluated.  In  addi¬ 
tion,  the  influence  of  the  crossing-trajectories  effect,  which  causes  ver¬ 
tical  and  horizontal  diffusivities  to  be  unequal,  is  taken  into  account. 
Throughout  our  analysis,  the  oil  is  assumed  to  be  in  the  form  of  small, 
noninteracting  particles:  the  important  question  of  the  process  by  which 
these  particles  are  formed  is  not  considered. 

Section  6  is  devoted  to  the  construction  of  a  suitable  computational 
method  to  determine  oil  concentrations  from  the  model  devised  in  §5. 
Since  the  diffusivities  are  very  small,  the  time  required  to  achieve  a 
steady  distribution  of  oil  is  found  to  be  too  large  and  the  spatial 
scales  too  small  to  allow  a  direct  finite  difference  simulation  to  resolve 
the  entire  region.  Instead,  we  resort  to  a  combination  of  analytical  and 
numerical  methods.  A  byproduct  of  this  method  of  solution  is  a  much 
clearer  picture  of  the  injection  and  loss  of  oil  to  SRZ's. 

Section  7  describes  the  results  of  our  computations  of  subsurface 
oil  dispersion. 

In  §8  an  attempt  is  made  to  describe  the  effects  of  LC  on  a  coherent 
layer  of  oil  floating  on  the  surface.  It  is  a  departure  from  the  main 
theme  of  the  work:  subsurface  motions  are  not  considered,  and  the  form 
that  the  oil  takes  is  not  that  assumed  in  the  rest  of  the  work.  A  simple 
model  is  developed  that  predicts  a  collection  of  oil  into  windrows,  and 
allows  an  estimate  of  the  amount  of  oil  present. 
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Finally,  in  $9,  we  discuss  the  conclusions  that  we  arrive  at,  and  we 
suggest  ways  that  the  information  we  have  found  may  to  be  put  to  use. 

3 . 1  Parameter  Ranges  for  Langmuir  Circulations  Computations 

The  complete  mathematical  model  for  Langmuir  circulations  will  be 

presented  in  $4;  it  requires  as  input  the  Stokes  drift  of  the  surface 
waves,  or  characteristic  wave  frequency,  the  wind  stress  applied  at  the 
ocean  surface,  and  the  eddy  viscosity  representing  the  effects  of  turbu¬ 
lence.  In  our  model,  we  represent  wind  generated  wave  fields  by  a  single 
representative  wavetrain  obtained  by  considering  known  properties  of 
empirically  determined  ocean  wave  spectra.  Since  the  spectra  depend  upon 
wind  speed  and  upon  the  wind  fetch  and/or  duration,  these  factors  will  be 
accounted  for  in  the  specification  of  the  wave  characteristics. 

The  wind  stress  will  also  be  related  to  wind  speed,  and  so  the  wind 
speed  Va,  and  the  sea  state  will  be  the  primary  environmental  inputs. 
Calculations  are  done  here  with  wind  speeds  Va  of  5,  10,  and  15ms-1. 

Three  choices  of  "ocean  mixed  layer  depth"  which  we  denote  by  H 
(taken  to  be  5,  10,  and  20  meters)  are  made  to  fix  the  maximum  depth  of 

penetration  of  LC  to  be  considered.  The  actual  mixed  layer  depth  may 
exceed  these  values,  and  in  fact  the  real  mixed  layer  plays  no  role  in  the 
calculation  except  that  we  expect  that  the  LC  depth  cannot  exceed  it. 
Thus,  for  each  value  of  H  and  wind  speed,  we  calculate  three  LC  fields, 
one  penetrating  to  depth  H,  and  two  others  penetrating  to  intermediate 
depths.  Altogether,  twenty-seven  LC  fields  are  calculated  in  this  way. 

We  now  describe  the  rationale  used  to  relate  each  of  the  input 
parameters  to  the  wind  speed. 

Wind  Stress 

The  water  friction  velocity,  u^,  is  a  measure  of  the  applied 
,  ,  2 

wind  stress,  ta,  which  is  taken  to  be  Ta=pwu*  ,  where  pw  is  the 
water  density.  Using  the  standard  estimates  (such  as  that  in  Leibovich 
(1977,  p.  740),  we  correlate  u*  with  wind  speed  through 

u*  3  Va/660  (3-1) 


27 


where  Va  is  the  wind  speed  a  few  meters  (5  or  more;  the  exact  height  is 
not  critical)  above  the  surface. 

Correlations  for  a  Fully-Developed  Sea 

The  spectrum  of  wind  generated  waves  is  sharply  peaked,  and  this 
spectral  peak  is  near  the  frequency  corresponding  to  the  waves  of  "signi¬ 
ficant  height"  (i.e.,  the  mean  of  the  1/3  highest  waves).  As  time  (wave 
duration  T,  in  cases  where  the  wind  is  switched  on  at  a  definite  time)  or 
distance  (wave  fetch  F,  in  steady  cases  with  wind  blowing  over  a  finite 
distance)  increases,  the  spectral  peak  moves  to  lower  and  lower  frequen¬ 
cies.  An  idealized  asymptotic  state  (as  duration  T  or  fetch  F  ♦  <*>)  in 
which  no  further  wave  growth  occurs  and  known  as  a  "fully-developed  sea", 
has  been  found  to  be  a  useful  hypothesis.  In  this  state  (the  existence  of 
which,  as  argued  by  Phillips  (1977),  is  conjectural),  the  wave  character¬ 
istics  (height  Hi/3,  frequency  o,  wavenumber  k,  phase  speed  c)  for  deep 
water  wind-driven  waves  may  be  related  to  the  wind  speed  as  follows: 


H] /3  3  aVa2/g 

(3-2) 

c  3  B  Va 

(3-3) 

and  we  assume  a  and  ic  are  related  to  c  through  the  linear  deep  water 
dispersion  relation 


c  *  o/tc  -  g/a  (3-4) 

where  g  is  the  acceleration  of  -ravity.  The  values  for  the  dimensionless 
constants  a  and  B  given  in  the  literature  vary.  Longuet-Higgins  (1969) 
takes  B  3  1,  a  *  0.25,  Stewart  (1967)  and  Neumann  &  Pierson  (1966)  also 
suggest  B  3  1,  but  take  a  3  0.2,  while  Bretschneider  (1966)  takes  B  3 
1.95,  a  3  0.283. 
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A  truly  Cully-developed  sea  may  never  be  approached  for  high  speed 
winds,  since  a  fetch  greater  than  the  circumference  of  the  earth  may  be 
required  (Bretschneider  (1966)),  and  this  must  be  recognized  when  esti¬ 
mates  using  the  above  values  for  a  and  g  are  used.  Real  wind  waves,  with 
finite  duration  T  and/or  fetch  F,  can  be  related  to  a  steady  imposed  wind 
field  Va  via  (3-3)  and  (3-4),  provided  a  and  g  are  taken  to  be  functions 
of  T  and/or  F.  For  a  developing  sea,  a  and  g  will  both  be  less  than  their 
fully  developed  values. 

In  the  LC  model  to  be  used,  the  natural  length  is  ic"1  ,  which  is 
A/2ir,  where  X  is  the  length  of  the  dominant  surface  waves.  A  relation 
between  tc  and  Va  is  found  from  (3-3)  and  (3-4), 

<=  B~2gVa"2  (3-5) 


Eddy  Viscosity 

The  mathematical  model  of  Langmuir  circulations  that  we  use  repre¬ 
sents  the  dynamical  effects  of  turbulence  on  the  coherent  convective 
motions  by  use  of  an  eddy  viscosity,  v-p.  This  is  the  input  parameter 
about  which  least  is  known,  and  it  is,  consequently,  the  most  difficult  to 
estimate.  Semi-empirical  estimates  which  correlate  eddy-viscosity  with 
surface  wave  characteristics  (see,  for  example,  Ichiye  (1967))  relate  vT 
to  the  significant  height  Hp/3  (average  height  of  the  highest-third  of 
the  waves  in  the  spectrum)  and  frequency  0  in  the  form 


v 


T 


CTH5/3CT 


(3-6) 


where  C-p  is  a  dimensionless  constant.  An  'Jstimate  of  this  form  is  inde¬ 
pendently  arrived  at  in  Leibovich  &  Radhakrishnan  (1977)  (equation  2)  by 

matching  theoretical  results  on  Langmuir  circulation  with  observat ions . 
.  .  3 

The  constant  C-p  in  that  case  is  0.7x10”  .  Two  other  estimates  are 
found,  also  independently,  in  that  paper;  using  correlations  appropriate 
for  a  fully-developed  sea,  Op  is  found  to  be  0.58x10"  and  1.5x10"  and 
the  mean  of  the  three  values  is  approximately  0.9xl0"3. 
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There  is  (inevitably,  due  to  the  scarcity  of  data,  and  the  intrinsic 
problems  of  using  and  of  specifying  an  eddy  viscosity  for  a  turbulent 
flow)  a  considerable  scatter  in  the  value  of  Cj,  although  the  order  of 
magnitude  is  consistent  in  the  estimates  cited  above,  and  is  consistent 
with  estimates  of  Cf  given  by  other  sources.  We  adopt  the  form  (3-6), 
and  in  the  section  on  "Langmuir  number",  we  choose  Cf  *  1.5x10“  .  It  is 
argued  there  that  the  final  results  in  the  Langmuir  circulation  problem 
ought  not  to  be  overly  sensitive  to  the  precise  value  of  Op  selected. 

The  Langmuir  Number 

In  water  of  constant  density,  the  LC  problem  can  be  reduced  to  a 
dimers  ionless  form  involving  one  dimensionless  parameter,  the  Langmuir 
number.  La.  This  parameter  depends  upon  sea  state,  eddy  viscosity,  and 
friction  velocity  as  follows: 

La  =  2(oovT^^/8h^3U*  (3-7) 

where  deep  water  gravity  waves  have  been  assumed,  and  where  the  character¬ 
istic  wave  frequency  and  height  are  taken  to  be  oq,  the  frequency  of  the 
spectral  peak,  and  Hp/3. 

Introducing  the  correlations  (3-1)  to  (3-4), 

La  =  1.32xl03CT3/2a2e”3  (3_g) 

which  is  independent  of  wind  speed.  If  we  adopt,  for  purposes  of  estima¬ 
tion,  a  *  0.25,  6  ■  1  (following  Longuet-Higgins  (1969)),  and  C-p  = 
1.5X10“  (representing  the  upper  value  of  those  discussed  below  (3-6),  we 
find 


La  «  0.005 


(3-9) 
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We  believe  that  this  is  a  reasonable  value  to  take  for  La,  despite  the 

uncertainities  in  the  parameters  Op,  a,  and  8.  Thus  we  adopt  (3-9)  in 
our  calculations. 

There  is  reason  to  believe  that  the  particular  form  of  the  Langmuir 
cells  may  not  be  sensitive  to  the  value  for  Vf  (hence  Oj)  in  any 

event.  Since  reasonable  estimates  generally  indicate  that  La  is  small, 
the  dynamics  are  nearly  inviscid.  The  eddy  viscosity  v-j>  is  likely  to 
play  a  significant  role  only  near  the  air-sea  interface,  where  it  serves 
as  the  mechanism  for  momentum  transfer  from  wind  to  water.  The  rate  of 
momentum  transfer  is  fixed  by  u*  only,  and  does  not  involve  v-j.  Below 
the  thin  surface  boundary  layer  (the  layer  in  which  v-j  is  dynamically 
significant),  momentum  transfer  is  accomplished  mainly  by  convective 
motions,  once  Langmuir  circulation  is  initiated.  Thus,  Vf  is  expected 

to  affect  the  onset  to  LC  convection,  but  once  convection  is  effectively 
initiated,  vT  should  primarily  determine  the  thickness  of  the  surface 

boundary  layer. 

LC  Length  Scale  Parameter 

The  inverse  of  surface  wavenumber,  or  x-1  ,  serves  as  the  length 
scale  in  the  LC  model.  Since  we  have  denoted  the  maximum  LC  penetration 
depth  by  H,  the  maximum  dimensionless  depth  allowed  is  xH  ,  For  a  given 
wind  speed  and  fetch,  one  may  compute  x;  together  with  H,  this  sets  a 
minimum  numerical  value,  xH,  of  the  (dimensionless)  computational  domain. 
The  depth  of  we  11 -developed  LC  corresponds  to  xH  on  the  order  of  unity. 
Very  little  qualitative  change  in  the  structure  of  LC  is  evident  beyond  xH 
*  3  to  4  in  the  published  simulations  of  LC.  Thus  we  have  arbitrarily  set 
the  maximum  dimensionless  value  of  the  depth  of  the  LC  field  calculations 
to  be  5.6.  This  is  large  enough  to  capture  all  known  Langmuir  features. 

We  fix  the  wavelength  of  the  surface  waves  under  consideration  by 
prescribing  xH  *  5:  in  this  way,  the  specification  of  H  fixes  x  for  the 
purposes  of  our  calculations  .  Given  the  wind  speed,  the  wave  field  under 
consideration  will  correspond  to  a  definite  fetch.  If  one  adopts  the 
empirical  spectral  development  model  to  be  presented  below,  this  fetch  can 
be  calculated.  (This  will  be  illustrated,  but  the  information  on  fetch  is 
not  essential  for  the  subsequent  work.)  More  important  is  the  fact  that 
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the  selection  of  x  also  determines  as  shown  below,  and  therefore 

the  sea  state. 

Development  of  the  Wind  Wave  Spectra  and  Stokes  Drift 

While  a  nearly  fully  developed  sea  may  be  established  in  light 
winds,  the  idealization  is  less  realistic  in  moderate  winds,  and  is  un¬ 
likely  ever  to  occur,  as  previously  mentioned,  in  heavy  winds.  If  we 
insist  upon  LC  penetration  depths  of  5,  10  and  20  meters  -  the  cases  dealt 
with  here  -  a  developing  sea  is  more  appropriate  than  a  fully-developed 
one,  except  at  the  lowest  wind  speed  (5m  s“*)  to  be  considered.  Further¬ 
more,  there  is  some  question  about  the  applicability  of  a  "fully-developed 
sea"  at  wind  speeds  as  high  as  15m  s-1 .  For  example,for  this  wind  speed, 
Bretschneider  (1966)  estimates  that  a  minimum  wind  duration  of  307  hours 
is  required  to  raise  a  fully-developed  sea.  We  note  that  even  a  very 
young  sea  will  generally  produce  LC  in  a  matter  of  a  few  tens  of  minutes. 
In  addition,  for  the  same  wind  speed  and  (5  *  1,  ( 3-5)predicts  x  “  0.04m, 
(length,  X,  of  surface  equals  144  meters),  so  the  length  scale  (x-1)  of  LC 
is  23m.  A  mixed  layer  with  a  depth  of  5m  is  thin  compared  to  the  depth  of 
LC  generation,  and,  once  LC  activity  is  initiated,  strong  convective 
activity  will  reach  to  depths  comparable  to,  Thus,  for  a  high  wind 

speed,  a  fully-developed  sea  predicts  wave  activity  that  is  too  vigorous 
and  unlikely  to  represent  sea  states  actually  realized. 

From  the  extensive  data  summarized  by  Phillips  (1977,  pp.  161-3), 
the  frequency  of  the  spectral  peak,  o0,  is  given  in  terms  of  the  fetch  F 
and  the  friction  velocity  u^a  of  the  airflow  above  the  water  by 

a0“2.2(g/u2*a)(gFu2*a)~1^  (3-11) 


and  the  mean-square  displacement  of  the  water  surface  is 


T2  -  1.6xlO-4(uZ*aF/g). 


(3-12) 
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If  the  development  of  the  spectrum  is  limited  by  wind  duration  T  rather 
than  fetch  F,  it  is  reasonable  to  interpret  (3-11)  and  (3-12)  by  replacing 
F  by  CgT,  where  Cg  is  the  group  velocity  of  waves  at  the  frequency  erg, 
or  cg=  l/2(g/o0). 

The  Stokes  drift  (see  S4)  for  the  ocean  contains  contributions  from 
all  the  waves  in  the  spectrum.  Our  model  replaces  these  by  a  single 
representative  wavetrain  with  frequency  ao  and  the  amplitude  of  the  "waves 
of  significant  height".  This  replacement  can  be  justified  by  considering 
the  effect  of  a  spectral  spread  on  the  drift,  following  (for  example)  the 
work  of  Kenyon  (1969). 

The  amplitudes  (half  height)  of  the  waves  of  significant  height, 
aj/3,  may  be  related  to  the  rms  amplitude  (3-12)  by  (Longuet-Higgins , 
1952) 


ai/3  =  1.42/S2  =  1/2Hj/3> 


(3-13) 


The  Stokes  drift  for  the  replacement  wave  is  (Phillips,  1977) 

2  2 

U8  ■  a  1/3k  ocoexp(2x0z) 


(3-14) 


where  cq  ■  a0/te0  i3  the  phase  speed  of  a  wave  having  frequency  o0  and 
wavenumber  Kq,  and  z  is  a  coordinate  measured  vertically  upwards  from  the 
mean  free  water  surface:  al/3ie0  i®  a  measure  of  the  slope  of  the 
replacement  wave.  Since  <o  i®  constructed  from  ao  by  the  deepwater 
dispersion  relation 


\ 


0O2  “  8*0 


(3-15) 


(3-10),  (3-11),  (3-12)  imply  a  wave-slope 


a  1  / 3*c o  *  0.087 


(3-16) 


which  is  independent  of  fetch  or  duration. 
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We  now  have  enough  information  to  assemble  a  table  summarizing  the 
set  of  LC  fields  generated.  Given  H  and  Va,  we  determine  Kq  from  (3-16) 
(in  considerations  of  the  effect  of  a  spectral  spread,  we  arrived  at  .088 
for  the  slope  parameter  rather  than  0.087,  and  it  is  the  former  value  we 
actually  have  used),  u*  from  (3-1),  and  La  from  (3-9).  From  these 
values,  the  convective  velocity  scale,  V,  of  the  LC  model  (see  S4)  can  be 
calculated.  This  is  given  by  (see  §4), 

V  *  ai/3u^(a/vT)^/^, 

Using  (3-6),  we  see  that  V  may  be  expressed  as 

V  =  u*/2/Ct  =*  12. 9u*  (3-17) 

if  we  use  the  value  suggested  for  Cf.  Thus  the  convective  velocity 
varies  linearly  with  wind  stress  if  the  parameterization  (3-6)  is  adopted; 
Thus,  using  (3-1),  we  relate  V  to  the  wind  speed  Va  through 

V  -  1.96*10-2Va.  (3~18) 

The  values  of  V  for  the  three  wind  speeds  in  our  calculations  are 
given  in  Table  3-1 


Table  3-1  -Entries  give  the  LC  convection  velocity  scale, 
V,  in  cm/s. 


Wind 

Speed  (in  meters/second) 

5 

10 

15 

(Case  A) 

(Case  B) 

(Case  C) 

1 

V 

1 

9.8 

19.6 

29.4 

Cases  1,  2,  and  3  correspond  to  depths  H  of  5,  10,  and  20  meters. 
The  length  scaling  factors  (kq_1)  and  wave  heights  H1/3  (-  2ax/3>  are 
the  same  for  all  wind  speeds,  and  are  given  in  Table  3-2. 
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Table  3-2.  -  LC  length  scaling  factors  Kg-*  in  meters, 
and  the  significant  wave  heights  for  cases  computed. 


Case 

K0-1  (in  meters) 

Hj/3  (in  :ra) 

1 

1 

17.6 

2 

2 

35.2 

3 

4 

70.4 

For  each  of  the  9  combinations  of  wind  speeds  and  maximum  depths  H 
of  Table  3.1,  we  have  computed  three  LC  fields.  As  described  in  §4,  the 
LC  model  predicts  a  continually  deepening  penetration  of  Langmuir  circu¬ 
lations.  For  each  of  the  nine  input  combinations,  w-»  select  LC  fields  at 
three  time  stages  in  this  continuing  development;  one  corresponds  to  the 
time  at  which  the  LC  have  penetrated  to  depth  H,  and  two  others  represent¬ 
ing  intermediate  stages  of  development.  So  far  as  the  dispersion  calcula¬ 
tion  is  concerned,  we  regard  each  of  these  fields  as  a  possible  steady 
realization  of  a  LC  field. 

It  is  worth  noting  that  the  cases  illustrated  here  for  the  highest 
wind  speed  and  lowest  H  values  correspond  to  very  immature  seas.  For 
example,  if  nterpreted  in  terms  of  wind  duration,  the  case  of  Va  =  15ra 
s'1 corresponds ,  according  to  the  estimates  presented  earlier,  to  a  dura¬ 
tion  of  only  5.25  minutes.  In  such  a  case,  our  dispersion  calculation  will 
probably  be  of  little  value,  since  the  seas  will  rapidly  become  much  more 
vigorous  than  assumed;  in  turn  this  means  that  the  LC  fields  can  only  be 
confined  to  the  assumed  small  depth  for  a  short  time,  and  soon  will  pene¬ 
trate  to  substantially  greater  depths.  We  could  (although  we  do  not  here) 
interpret  our  dimensionless  LC  fields  in  a  different  way  that  avoids  this 
loss  of  significance.  We  could  fix  <g  by  relating  it  to  fetch  via  (3-11) 
and  (3-15)  for  example,  or  by  arbitrarily  fixing  Hj/3  and  working  back¬ 
wards  to  *0  using  (3-16). 
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4.  MATHEMATICAL  MODELS  FOR  THE  SIMULATION  OF  LANGMUIR  CIRCULATION, 

STOMMEL  RETENTION  ZONES,  AND  OIL  DISPERSION 

The  wave-current  interaction  theories  of  Langmuir  circulations 
describe  the  water  motions  by  the  averaged  Navier-Stokes  equations. 
Averages  are  taken  over  times  comparable  to  the  period  of  surface  waves, 
and  so  the  motions  described  by  the  resulting  set  of  equations  have  sur¬ 
face  wave  frequencies  filtered  out.  The  required  averaging  process, 
(which  is  not  straightforward),  is  given  in  Craik  and  Leibovich  (19'S)  and 
more  systematically  in  Leibovich  (1977a).  The  theory  can  also  be  derived 
from  the  generalized  Langrangian  mean  formalism,  an  exact  theory  of  wave- 
current  mean  flow  interactions  due  to  Andrews  and  McIntyre  (1978):  a  deri¬ 
vation  of  the  LC  governing  equations  by  this  route  is  given  by  Leibovich 
(1980).  In  the  theory,  it  is  assumed  that  incoherent  (turbulent)  motions 
may  be  described  by  a  constant  eddy  viscosity. 

Following  the  presentation  of  the  LC  model  a  brief  discussion  of  its 
numerical  solution  in  sections  4.1  and  4.2,  we  turn  to  a  description  of 
Stommel  retention  zones  in  §4.3;  the  diffusion  model  to  be  used  for  the 
dispersion  calculations  is  then  is  then  presented  in  §4.4.  The  important 
question  of  the  justification  for  the  use  of  a  diffusion  model  for  the 
dispersion  is  given  an  entire  section  (5). 


4.1  Governing  Equations  for  Langmuir  Circulations 


The  equations  governing  the  formation  of  Langmuir  circulations  in 
water  of  constant  density  are,  according  to  Leibovich  (1977a), 


v,t  +  (  v*V  )va‘-Vx  +  u8  x(curl  v  )  +  V  2 


V»v  -  0 


(4-2) 


where  a  comma  denotes  partial  differentiation,  v  is  the  Eulerian  mean 
velocity  vector  in  the  water,  it  is  a  scalar  modified  pressure,  and  u8 
is  the  Stokes  drift  due  to  the  water  waves.  In  the  absence  of  u8  ,  and 
if  Vj  is  interpreted  not  as  an  eddy  viscosity  (as  we  shall  do)  but  as 
the  molecular  velocity,  then  (4-1,2)  are  the  Navier-Stokes  equations  for 
the  conservation  of  momentum  and  mass  of  an  incompressible  fluid.  The 
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rectified  effects  of  the  surface  waves  appear  only  in  the  form  of  the 
"vortex  force" 

u8  * (curl  v  ).  (4-3) 


Here  the  Stokes  drift  of  surface  waves  of  small  slope  is  given  by 
(Phillips,  1977) 


Uj,  dt*V) 


(4-4) 


where  u„  is  the  wave  velocity  vector  and  the  overbar  represents  a  time 
average  over  intervals  extending  over  several  surface  wave  periods.  In 
this  work,  we  assume  that  us  is  in  the  direction  of  the  wind  stress, 
taken  to  be  the  x-direction,  and  that  us  varies  only  with  depth  z 
(measured  vertically  upwards  from  the  mean  free  surface).  The  water  depth 
is  taken  to  he  infinite,  and  so  the  ocean  lies  in  -»<z<0.  The  assumptions 
made  about  u8  are  consistent  with  a  random,  wind-generated  sea  over 
restricted  fetch  intervals.  We  make  the  further  assumption  that  we  may 
replace  u8  by  the  Stokes  drift  of  a  "replacement  wave"  (see  §3),  and 
take  (3-14)  for  |  u8  j . 

The  ocean  is  assumed  to  be  at  rest  initially,  and  to  begin  to  move 

2 

in  response  to  an  applied  wind  stress  t  =  pu*  and  surface  wave  field 
imposed  at  time  t  =  0.  If  we  take  (u,v,w)  to  be  the  (x,  y,  z)  velocity 
compoments,  then  the  boundary  and  initial  conditions  on  the  velocity 
vector  v  (x,y,z,t)  are 

U,z=  U  */vT,  U,z  =  w  a  0  on  z  =  0  (a) 


(u,v,w)  +  0  as  z  -  °° 


(b)  (4-5) 


(u,v,w)  =  0  at  t  =  0  (c) 

The  problem  is  made  dimensionless  by  the  following  scaling  scheme: 
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Length  scale:  Kg-  (a) 

velocity  scale  in  x  direction  =  U  ■  u*  /vx<g  (b) 

velocity  scale  in  the  y  and  z  directions  (convective  scale) 

=  V  =  a\/ 3U*(oo/\*t^^  (c) 

time  scale:  (<oV)-1  (d)  (4-6) 

2  2 

modified  pressure:  u  a  .,(ag/v  )  ,  „ 

l/J  i  (e) 

Stokes  drift:  a^,_tc0o0. 

1/3  0  0  (f) 


When  substituted  into  (4-1),  (4-2)  and  (4-5),  the  dimensionless  problem  is 
(we  use  the  same  symbols  for  dimensionless  variables  and  their  dimensional 
counterparts  since  no  confusion  is  likely) 


r,  t  +  (  )v=u8Vu-Vx+La  V ■w  (a) 

7-v  *  0  (b) 


(4-7) 


2 

and  (4-5)  remains  the  same  except  that  u,2mu*  /vx  is  replaced  by 
u,z  3  1  in  (4-5a) .  The  parameter  La  is  the  Langmuir  number  defined  by 


La  *  kqvt/V. 


(4-8) 


The  solution  to  the  problem  as  posed  is  inherently  time-dependent, 
since  the  force  applied  at  the  water  surface  is  resisted  only  by  the  fluid 
inertia.  A  steady  motion  is  obtainable  only  by  taking  into  account 
effects,  such  as  the  Coriolis  acceleration  due  to  the  Earth's  rotation 
(important  for  problems  where  the  motion  persists  for  times  larger  than 
about  one  quarter  pendulum  day),  a  finite  bottom  (which,  under  most 
circumstances  is  probably  irrelevant),  or  some  less  obvious  sink  for 
momentum.  We  accept  the  present  problem,  and  calculate  the  evolving 
mot  ion. 

It  is  worth  noting  that  a  purely  rectilinear  solution  is  possible 
(Leibovich,  1977a)  in  which  v  *  w  ■  0  and  u  varies  only  with  z  and  t. 
This  motion  is  unstable  for  La  <0.66  (Leibovich  and  Paolucci  (1981)), 


\ 
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however,  and  the  expected  form  of  the  instability,  according  to  the  theo¬ 
retical  work  of  Leibovich  and  Palucci  (1980b),  is  rolls  (disturbances 
independent  of  x)  ,  consistent  with  the  long  windrow  patterns  actually 
observed  in  the  ocean.  To  get  such  an  instability  numerically,  one  need 
only  calculate  the  rectilinear  solution  out  to  some  arbitrary  time  level, 
and  then  superimpose  at  that  time  a  very  small  perturbation-such  as  a  weak 
spanwise  sinuosity-and  continue  the  computation.  If  La  <0.66,  the  convec¬ 
tive  rolls  will  grow  with  time  to  produce  well  formed  simulated  Langmuir 
circulations.  The  next  section  discribes  the  numerical  method  used. 

4 . 2  Numerical  Method 

The  solutions  were  carried  out  by  the  explicit  finite  difference 
algorithm  developed  by  Radhakr ishnan  (1975)  and  modified  by  Paolucci 
(1979).  A  description  of  the  method,  which  is  applied  to  the  x-indepen- 
dent  version  of  (4-7)  in  vorticity-strearafunction  form,  is  given  by  Leibo¬ 
vich  and  Paolucci  (1980a),  and  will  not  be  reproduced  here. 

It  is  necessary  to  choose  a  finite  spanwise  dimensionless  width  for 
the  computational  domain,  and  our  computations  were  done  for  0  y  4. 
This  is  sufficiently  large  to  encompass  Langmuir  circulations  with  a 
windrow-to-windrow  dimensional  distance  of  8<o_1  or  smaller.  For  the 
cases  computed  this  allows  the  description  of  LC  with  windrow  spacings  up 
to  32  meters. 

The  spatial  mesh  used  was  Ay  =  Az  =  0.1,  which  implies  a  cell  Rey¬ 
nolds  number  of  around  10.  The  discussion  of  numerical  diffusion  given  by 
Leibovich  and  Paolucci  (1980)  applies  here  as  well. 

By  placing  the  problem  in  dimensionless  form,  we  are  able  to  reduce 
all  cases  to  be  run  (see  Table  3-1)  to  a  single  numerical  problem.  The 
initial  value  problem  is  calculated  until  the  convective  activity  has 
penetrated  to  a  dimensionless  depth  of  5.6  (see  §3.1)  and  complete  sets  of 
the  three  velocity  components  are  stored  at  numerous  intermediate  time 
levels.  From  this  collection,  three  data  sets  are  selected,  the  final  one 
and  two  representative  sets  at  intermediate  times.  Each  of  these  is 
converted  to  nine  dimensional  velocity  fields  by  use  of  the  scaling 


factors  (4-6)  and  Tables  3-1  and  3-2.  This  produces  27  different  velocity 
fields,  each  of  which  is  regarded  as  a  possible  realization  of  a  Langmuir 
circulation  convective  motion.  The  three  basic  dimensionless  LC  fields 
are  shown  in  figures  4-1,  2,  3. 

As  an  example  of  a  dimensional  realization,  we  note  that  for  case  C, 
with  a  windspeed  10m  s-*,  figure  4-3  describes  a  case  for  which  the  dis¬ 
tance  between  windrows  is  16  meters.  Each  of  these  is  stored  on  magnetic 
tape,  and  converted  to  a  Stommel  retention  zone  as  described  in  the  next 
subsection. 

4 . 3  Calculation  of  Stommel  Retention  Zones 

The  vertical  velocity  in  LC  vanishes  at  the  surface  ,  and  grows  in 
magnitude  with  depth  below  convergence  lines.  If  a  buoyant  particle  is 
introduced  at  a  point  below  the  convergence  line  where  downwelling  speeds 
exceed  the  (positive)  buoyant  terminal  velocity  of  rise,  VT,  the  parti¬ 
cle  will  be  carried  down.  On  the  other  hand,  the  particle  will  eventually 
be  swept  laterally  away  from  its  initial  position  into  a  region  where 
water  downwelling  is  less  than  Vf,  and  the  particle  will  tend  to  rise 
back  towards  the  surface.  A  fascinating  result  of  Stommel's  (1949)  (in  a 
paper  in  which  an  artificial  velocity  field  was  used  to  show  that  LC  can 
prevent  heavy  particles  from  sinking)  is  that  a  buoyant  particle  can  be 
trapped  forever  in  a  zone  some  distance  below  the  surface  if  it  is  inject¬ 
ed  into  that  zone  and  if  turbulence  is  ignored.  A  subsurface  zone  in 
which  particles  can  become  trapped  will  be  called  a  "Stommel  retention 
zone"  or  SRZ  here  .  These  zones  will  play  a  major  role  in  the  oil  disper¬ 
sion  calculation  method  introduced  in  §4.4. 

We  continue  to  use  the  same  coordinate  system  of  and  suppose  an  LC 
field  has  developed.  Then  the  horizontal  and  vertical  velocity  compon¬ 
ents,  v  and  w  may  be  given  in  terms  of  a  streamfunction  in  tp,  or 

v  -  i|»,z,  w  -  -  <p,y.  (4_8) 

If  we  assume,  as  Stommel  did,  that  the  velocity  of  the  oil  relative  to  the 
surrounding  water  is  always  in  the  vertical  direction  with  constant  speed 
Vf,  then  the  trajectory  of  an  oil  particle  y  *  Y(t),  z  ■  Z(t)  is  given 
by 
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Y,t  -  v  =  *>z  (a) 

Z,t  »  v  +  VT  >  -  (tJ>-VTy),r  (b) 


(4-9) 


If  the  fluid  motion  is  steady  in  time,  then  3i|i/3t  ■  0,  and  trajec¬ 
tories  of  oil  particles  in  the  cross-plane  determined  by  (4 . 9)coinc ide 
with  level  curves  of  the  function 


¥(y,z)  »  4>(y , z )  -  VTy. 


(4-10) 


It  is  not  hard  to  see  that,  if  »  *  Vj  <  0  for  given  Vf  anywhere  on  a 
downwelling  plane,  then  a  Stommel  retention  zone  exists  for  all  particles 
with  terminal  velocity  of  rise  less  than  or  equal  to  V^.  These  zones 
are  nested  regions  for  which  the  contours  of  ¥  are  closed.  Each  of  these 
nested  closed  contours  (if  any  exist)  is  enclosed  within  its  largest  mem¬ 
ber,  which  forms  the  boundary  of  the  SRZ.  The  disposition  and  extent  of 
an  SRZ  depends  critically  upon  Vf,  the  LC  convective  velocity  scale  V 
and  length  scale  x0-1  (since  the  fluid  st rearafunc t ion  is  proportional  to 
V<0-').  If  C  is  the  boundary  contour  of  an  SRZ  for  a  given  value  of  Vj 
in  a  given  LC  field,  than  any  particle  with  this  value  of  Vj  or  smaller 
placed  inside  C  will,  according  to  this  model,  remain  trapped  within  C 
forever.  If  VT  is  reduced,  the  SRZ  is  enlarged;  if  increased  the  SRZ  is 
reduced. 

The  level  lines  of  ¥  are  plotted  in  Figure  4-4  for  Vf  ■  0.84  cra/s, 
Va  =  10  m/s,  Hj/3  «  0.7  m.  Note  that  the  SRZ  is  very  extensive,  and 
distribute  oil  to  depths  exceeding  20  meters.  Retention  zones  for  a 
number  of  terminal  velocities  are  given  in  the  Appendix:  the  examples 
cover  all  twenty-seven  LC  fields  that  we  have  computed. 
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4.4  The  Diffusion  Model  for  Oil  Dispersion 

We  attempt  to  describe  the  concentration  of  oil,  that  is,  the  volume 
fraction  of  the  water  volume  occupied  by  oil  at  any  position  and  time; 
this  is  equivalent  (Tennekes  and  Lutnley,  1972);  Monin  and  Yaglora,  1971)  to 
asking  for  the  probability  that  a  single  oil  particle  may  be  found  in  the 
subvoturae  centered  on  the  position  of  interest.  The  concent  rat  ion,  C, 
will  be  assumed  to  be  controlled  by  the  following  diffusion  equation 

c,t  ♦  v-V  C  +  VTC,s  =  (DjC.yKy  ♦  (DlC>2),2  (4-11) 

where  we  take 

C  -  C(y,z,t)  (4-12) 

▼  is  the  velocity  vector  of  the  mean  motion  (here  taken  to  be  Langmuir 

circulations),  and  Vj  is  the  terminal  velocity,  positive  for  positively 

buoyant  oil  particles.  The  parameters  Dj  and  D2  are  di f fusivi t ies  result¬ 
ing  from  turbulent  buffeting  of  the  oil  and  are  in  general  different  for 
the  horizontal  (D^ )  and  the  vertical  (D2)  directions;  they  may  depend  upon 
depth. 

The  derivation  of  this  equation  is  discussed  in  Luraley  (1978b)  and 
the  validity  of  its  use  for  oil  particles  in  water  is  considered  in  §5. 

In  calculations  that  we  will  do,  it  will  be  assumed  that  the  amount 
of  oil  in  the  region  of  the  ocean  considered,  is  fixed.  That  is,  no  new 
oil  is  being  injected,  and  our  job  is  to  determine  how  the  fixed  amount  of 

oil  is  distributed.  The  total  amount  of  oil  in  the  volume  V  of  the  ocean 

being  considered  is 


V  *  /c  dxdydz  . 

v  (4-13) 

Its  rate  of  change  with  time  is  zero,  by  hypothesis.  The  Langmuir  circu¬ 
lation  cells  described  in  sections  (4.1),  (4.2)  are  periodic  in  the 

cross-wind  horizontal  direction  y  and  are  independent  of  x.  Consider  a 
volume  V  of  unit  distance  in  the  x-direction,  and  consisting  of  one  period 
(0,Yi)  of  the  LC  in  the  spanwise  (y)  direction,  and  of  infinite  depth. 
Since  the  spanwise  distribution  of  oil  is  determined  strictly  by  the  LC 
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cell,  C  will  have  the  same  periodicity,  and  C,y  vanishes  of  the  bounda¬ 
ries  of  the  ceil,  say  at  y*0  and  y*Yi  .  Integrating  (4-11)  over  this 
volume  (0  z  >-°° ,  0  x  _<  1,  0  jCy^Yj)we  find 

V,t  -  /*l[D2(0)C,z(y,0,t)  -  VTC(y,0,t)]dy 


since  we  must  have  w(y,0,t)  =  0  in  the  Langmuir  circulations.  Since  V 
does  not  change  with  time,  the  integral  must  vanish.  The  integral  is  the 
flux  of  oil  into  the  ocean  across  the  mean  free  ocean  surface,  so  we  take 
as  one  boundary  condition 

VTC  -  D2  C , z  “  0  at  z  ■  0,  (4-14) 

and  we  also  have  the  additional  conditions 

C,y«  0  at  y  «  O.Yj  (4-15) 

C  *  0  (J  t  ‘ 

(4-16) 


These  boundary  cnditions  hold  for  unsteady  as  well  as  steady  conditions, 
although  we  will  concern  ourselves  only  with  the  time-independent  case. 

The  solution  of  the  problem  posed  by  (4-11),  (4-14),  (4-15)  and 
(4-16)  will  be  dealt  with  in  56,  Here  we  make  two  observations.  First, 
the  level  of  C  is  not  determined;  if  we  find  a  solution,  any  multiple  of 
it  is  also  a  solution.  The  magnitude  of  C  is  determined  by  the  additional 
condition  (4-13)  which  is  just  an  initial  condition  for  the  unsteady  prob¬ 
lem  whose  steady  limit  we  seek. 

Second,  the  role  of  the  Stommel  retention  zones  can  now  be  further 
clarified.  If  we  make  (4-11)  dimensionless  using  LC  velocity  scale  V,  and 
LC  length  scale  <q_1 ,  then  we  may  write  the  steady  version  of  (4-11)  in 
dimensionless  form  as 
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vC  > y+  (W  +  VT/V)  C,z  =  (e 1 C , y )y+  (e2c>z),2  (4_17) 

where 

G1  =  |C0D1/V>  e2  e  k0^2 ^ 

(4-18) 

and  otherwise  we  use  the  same  symbols  for  dimensional  and  dimensionless 
variables . 

We  find  in  55  that  Cj  and  £2  are  small  (of  order  10"3 ) .  Consequent¬ 
ly,  it  is  clear  that  the  reduced  equation  found  by  neglecting  C\  and  C2 
altogether  will  play  an  important  role.  Returning  now  to  dimensional 
variables,  the  reduced  equation  is 

v  C,y+  (w  +  VT)  C,z  »  0,  (4-19) 

thus,  as  usual  in  problems  of  boundary  layer  character,  the  leading 
approximation  satisfies  an  equation  of  reduced  order  and  we  are  not  able 
to  enforce  all  boundary  conditions. 

The  general  solution  to  (4-19)  is 

c  -  C«>  (4-20) 


where 


v  =  ? ,z>  w  ♦  vx=  -  T,y  (4-21) 

so  that  C  is  constant  on  lines  of  constant  ¥.  These  lines,  in  turn,  are 
the  level  curves  found  in  the  previous  section;  a  closed  ¥  *  constant  line 
defines  an  SRZ. 

The  restoration  of  the  diffusivities ,  and  the  enforcement  of  the 
boundary  conditions  is  left  to  56. 
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5.  THE  MATHEMATICAL  MODEL  FOR  DISPERSION  OF  OIL  IN  TURBULENT  FLOW 

Lumley  (1978)  examines  the  conditions  under  which  one  may  treat  the 
dispersion  of  small  particles  in  a  turbulent  fluid  by  the  use  of  a  diffus¬ 
ion  equation.  In  the  last  part  of  §4,  we  introduce  the  diffusion  model; 
the  limits  of  its  validity  are  examined  in  detail  in  the  present  section. 

In  Lumley' s  approach,  steady  state  d i f fus ivi t ies  are  determined,  assuming 
that  the  droplet  drag  law  is  linear.  The  principle  problem  is  to  deter¬ 
mine  whether  the  assumption  of  a  linear  drag  law  is  justified,  and  to 
estimate  the  terminal  velocity  corresponding  to  it.  Most  droplets  of 
interest  are  too  large  to  be  in  the  Stokes  regime,  so  that  a  linear  law 
results  only  from  linearization  about  a  mean  velocity.  This  results  in 
slightly  different  drag  coefficients  in  the  horizontal  and  vertical  direc¬ 
tions.  Below  we  will  consider  various  aspects  of  the  question  of  limita¬ 
tions  on  particle  size  which  will  permit  a  linear(ized)  drag  law  for 
particles  in  a  turbulent  fluid,  and  also  subject  to  wave  accelerations  and 
a  large-scale  mean  flow. 

5 . 1  Permissible  Fluctuating  Relative  Reynolds  Numbers 

We  may  write  the  instantaneous  drag  force  as  (U£/U)f(U),  where  U 
is  the  magnitude  of  the  total  instantaneous  relative  velocity,  and  is 
the  vector  value  of  the  total  instantaneous  relative  velocity.  This  pre¬ 
sumes  that  the  wake  is  steady,  and  that  changes  in  the  oncoming  flow  are 
sufficiently  slow  that  the  boundary  layers  and  wake  have  time  to  adjust  to 
the  changed  conditions.  It  implies  that  the  drag  is  always  aligned  with 
the  flow,  although  it  is  a  non-linear  function  of  the  magnitude.  This  is 
discussed  at  length  in  Lumley  (1978).  The  function  f  is  given  by  (to 
second  order) 

f  -  (U/as)[l  ♦  (3/ L6)R]  (5-1 ) 
where  ag  is  the  Stokes  value  of  the  time  constant,  given  by  Vfg/g, 
where  Vjg  is  the  Stokes  value  of  the  terminal  velocity  given  by  equation 
(5-6).  Now,  we  may  expand  the  expression  (U{/U)f(U)  in  a  Taylor  series 
about  the  value  *  (0,0, Vj)  keeping  terms  through  second  order. 
The  derivatives  are  given  in  general  by 

I  -  (f/U)(«ij  -  UiUj/U2)  ♦  f'UiUj/U2 
It  *  (f'/U  -  f/U2)(6(jUk/U  ♦  <5  ik.U  j/u  ♦  «jkUi/U  - 

3UiUjUk/U3)  +  f"UjUjUk/U3  (5-2) 
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where  we  have  indicated  by  I  and  II  the  first  and  second  derivatives 

respectively*  and  by  f'  and  fM  the  first  and  second  derivatives  of  f.  If 

2 

we  write  f  =  U/ag  +  BU  ,  then  in  a  horizontal  direction  (say,  the 
l-direction)  we  have  as  a  requirement  that  the  second  derivative  be  much 
smaller  than  the  first  derivative,  that 

“UjujB  <<  f/U  (5-3) 


substituting  (5-1),  this  translates  to 

( 3/ 16) R'  «  1  ♦  ( 3/ 16) R  ( 5-4) 

where  R’  is  the  fluctuating  relative  Reynolds  number.  If  we  take  R'  < 
0.218R,  then  (5-4)  will  be  satisfied,  with  the  left  side  being  not  above 
20%  of  the  right  side,  if  R  *  60,  the  largest  Reynolds  number  which  will 

permit  a  steady  wake.  In  the  vertical  direction  (the  3-direction),  if  we 

take  the  turbulence  approximately  isotropic,  we  obtain 

(6/16)R‘  «  1  +  (6/ 16)R  (5_5) 

Taking  R'  _<  0.21R  will  again  assure  that  the  first  neglected  terra  will  be 
no  more  than  20%  of  the  terra  retained.  With  R  *  60,  this  suggests  that  R' 
=*  12.5  is  a  safe  value. 

5 . 2  Limits  on  Particle  Size 

The  expressions  given  in  Lumley  (1978)  pp,  297,  301  must  be  modi¬ 
fied,  because  the  ratio  0o/pw  is  not  large  relative  to  unity,  as  was 
assumed  there.  The  modified  forms  are  (obtained  from  Lumley  1957): 

Vt  “  -d2g(po/Pw  ~  1  >/  18v  (5-6) 

R  -  -d3g(po/pu  -  l)/18v2  (5-7) 
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8  =  (2p0/P„  ♦  l)(d/n)3/38 


(5-8) 

R  »  -(p0/pw  -  1H2p0/pu  ♦  l)(d/n)5(2g/u2)/  1352Rtl/4  (5.9) 

where  the  notation  is  the  same  as  that  in  Luraley  (1978):  d  is  droplet 
diameter,  g  is  the  acceleration  of  gravity,  the  subscripts  0  and  w  refer 
to  oil  (droplet)  and  water  respectively,  v  is  the  kinematic  viscosity,  n 
i  is  the  Kolmogorov  microscale,  l  the  integral  scale  of  the  turbulence  and 

|  u  is  the  r.ra.s.  turbulent  fluctuating  velocity.  Equation  (5-6)  gives  the 

Stokes  terminal  velocity,  equation  (5-7)  the  relative  Reynolds  number 
based  on  the  terminal  velocity,  equation  (5-8)  the  relative  Reynolds  num¬ 
ber  (based  on  the  fluctuating  relative  velocity)  supposing  that  the  drop¬ 
let  is  approximately  Lagrangian  (i.e.  -  follows  the  fluctuating  velocity) 
;  in  its  behavior,  and  (5-9)  the  relative  Reynolds  number  (based  on  the 

•  fluctuating  relative  velocity)  supposing  that  the  behavior  of  the  particle 

j  is  approximately  Eulerian  (i.e.  -  rises  through  the  turbulence  being 

little  influenced  by  it).  Since,  for  a  given  size  oil  droplet,  we  do  not 
j  know  £  priori  what  its  behavior  is,  we  must  calculate  the  various  relative 

1  Reynolds  numbers  to  determine  its  behavior. 

’  For  a  typical  situation  we  take  a  mean  wind  in  the  atmosphere  of  10 

j  ra/s,  and  a  wave  height  of  2  ra.  Supposing  that  the  wind  stress  goes  into 

the  turbulence  primarily,  we  have 

u2  *  (pa/pw)v^2  *  (pa/pw)<v*/ua>2ua2  (5-10) 

where  the  subscripts  a  and  w  refer  to  air  and  water  respectively.  v+ 
is  the  friction  velocity  in  the  air.  If  we  take 

(vA/Ua)2  -  10-3,  (pa/pu)  »  10-3  (5-11) 
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then 


2 


2  ,  2 
ii  /  s  . 


We  take  1  to  be  of  the  order 
w 

which  is  a  reasonable 


u  “  10 

This  gives  Rg  *  2*10 

=  1 . I9x 10‘3  tn.  This  is  also  typical. 
Computing  values  of  V^/u  for  two  values  of 
left),  we  obtain  for  various  d  (across  the  top)) 


of  the  wave  height, 
value,  while  o  s 

the  rat io  Pq/pw  (°n 


100pm 

300pm 

lOOOura 

3000  pm 

0.70 

1 ,63x I0_1 

1.47 

1 ,63xlo' 

1  .47  xlO2 

0.95 

2.72x IO'2 

2 .45x  10-1 

2.72 

2.45x10’ 

If  the  terminal  velocity  19  large  relative  to  the  r.m.s.  fluctuating 
velocity,  we  may  regard  the  droplet  as  essentially  Eulerian,  and  vice  ver¬ 
sa.  We  see  that,  with  Po/Pw  =  0*7,  100  um  particles  and  smaller  may 
be  consideied  Lagrangian,  while  particles  of  1000  pm  and  larger  are 
Eulerian;  at  Pq/pw  =  0.95,  particles  of  300  pm  and  smaller  are 

Lagrangian,  while  those  3000  pm  and  larger  are  Eulerian. 

If  we  now  apply  the  Eulerian  and  Langrangian  criteria  for  the  rela¬ 
tive  Reynolds  number  to  be  below  12.5,  we  find  as  limiting  diameters 

Oo/pu  0.7  0.95 

Eulerian  1277  pm  1761  pm 

Lagrangian  6930  pm  6521  pm 

In  fact,  the  full  relationship  i9  instructive,  and  is  graphed  below: 
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100pm  300pm  1000pm  3000pm  10000pm  Ln(d) 

The  curve  marked  *,  having  a  slope  of  -2,  is  the  dividing  line  between  the 
Eulerian  and  Lagrangian  cases;  Lagrangian  cases  lie  to  the  left  of  the 
line,  and  Eulerian  to  the  right  of  the  line.  The  curve  designated  by  @  is 
the  Eulerian  criterion;  to  the  left,  the  fluctuating  relative  Reynolds 
number  is  12.5.  The  curve  designated  &  is  the  Lagrangian  criterion;  to 
the  left  the  fluctuating  relative  Reynolds  number  is  less  than  12.5. 

It  can  be  seen  that  the  Eulerian  and  Lagrangian  criteria  do  not 
cross  until  pg/pw  is  of  the  order  of  0.999,  and  that  both  cross  the 
Eulerian/Lagrangian  dividing  line  at  roughly  this  point,  at  a  particle 
diameter  of  order  4500  pm  at  this  density  ratio.  Hence,  for  our  densi¬ 
ties,  the  Eulerian  criterion  is  the  operative  one. 

So  far  as  the  Reynolds  number  based  on  the  terminal  velocity  is  con¬ 
cerned,  taking  the  limiting  Eulerian  diameter,  we  obtain  a  value  of  340for 
the  case  Pg/Pw  «  0.7,  and  a  value  of  149  for  the  case  p0/p  u  ■> 
0,95.  Since  these  values  are  well  into  the  non-linear  range,  the  drag 
will  be  higher,  and  the  terminal  velocity  lower,  than  the  Stokes  value; 
hence,  the  true  Reynolds  number  will  be  lower.  We  recall  that  Van  Dyke 
0975)  gives  R  ■  60  as  the  approximate  limit  for  steady  flow.  We  must 
examine  the  influence  of  non-linearity  on  the  Reynolds  number. 
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In  Van  Dyke  (  1975)  we  can  find  an  expression  for  the  drag  coeffi¬ 
cient  on  a  sphere: 


CD  =  (6*/R)[l  +  ( 3/8 ) R  +(9/40)R2lnR  +  0(R2)] 


(5-12) 


2  2 

where  Cp  =  D/pU  a  ,  where  D  is  the  drag  and  a  the  droplet  radius.  R  3 
Ua/v,  the  Reynolds  number.  The  terminal  velocity  is  given  by  equating  the 
buoyancy  and  the  drag.  If  we  agree  to  keep  only  the  first  two  terms  of 
the  drag  law,  and  if  we  revert  to  our  definition  of  R  based  on  droplet 
d i ame ter,  we  obtain 


RS  »  R[1  +  (3/ 16) R]  (5-13) 

where  Rs  is  the  Stokes  value.  If  we  take  R  =  60  as  the  limit  for  a 
steady  wake,  we  obtain  a  value  of  Rs  of  735.  Using  this  value,  we  ob¬ 
tain  limiting  diameters  of  d  *  1650  pm  for  Pq/pw  =  0.7  and  d  ■  3000  pm 
for  P q/ Pw  =  0.95. 

We  may  conclude  on  the  basis  of  this  paragraph  and  5.2  that  parti¬ 
cles  of  roughly  1500  pm  and  smaller  may  safely  be  considered  as  having 
steady  wakes,  and  linearization  around  the  terminal  velocity  is  perraiss- 
able.  1500  pra  is  also  approximately  the  limit  for  uniformity  of  the  flow 
field,  since  n  is  approximately  of  this  order. 

5 .4  Effect  of  Excessive  Particle  Size 

As  we  consider  larger  and  larger  droplets,  a  number  of  possibilities 
must  be  considered:  for  example,  in  computing  the  relative  velocity,  it 
was  assumed  that  the  frequency  corresponding  to  the  particle  response  time 
was  above  the  spectral  cutoff  of  the  turbulence.  If  this  is  no  longer 
true,  some  changes  must  be  made  in  the  various  estimates.  If  we  are  in 
the  non-linear  regime,  we  can  compute  the  horizontal  and  vertical  time 
constants  for  the  particle  motion,  and  we  find  that  in  the  worst  case,  the 
diameter  at  which  the  time  constant  would  be  the  same  as  the  spectral 
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cutoff  would  correspond,  for  a  density  ratio  of  0.7,  to  a  diameter  of  6190 
pro,  or,  for  a  density  ratio  of  0.95,  to  a  diameter  of  5200  urn.  Since  the 
maximum  size  of  particles  we  can  consider  (from  a  wake-stability  point  of 
view)  are  considerably  smaller  than  these,  we  may  compute  the  relative 
velocity  supposing  that  the  particle  characteristic  frequenvy  is  greater 
than  the  spectral  cutoff. 

Hence,  even  for  these  quite  large  particles  of  up  to  3000  um,  there 
is  no  cause  for  concern  that  any  of  our  assumptions  will  be  violated. 
Finally,  we  may  note,  however,  that  even  for  larger  particles,  that  might 
violate  our  assumptions,  there  will  be  no  ultimate  effect  on  the  diffusiv- 
ity.  In  the  crossing-t ra jector ies  formula  of  Csanady  (1963)  for  example, 
the  effect  of  too  targe  a  particle  time  scale  will  be  to  change  the  equi¬ 
valent  fluid  point  di f fusivi ty ,  which  will  be  that  resulting  from  a  trun¬ 
cated  spectrum.  However,  the  integral  in  question  depends  primarily  on 
the  energy-containing  range  of  wavenumbers,  and  is  relatively  little 
affected  by  the  high  wavenumber  end  of  the  spectrian  -  at  any  reasonable 
integration  time  (say,  for  times  greater  than  two  integral  scales)  the 
change  in  the  cut-off  will  have  essentially  no  influence.  Hence,  the 
particle  diffusivity  will  be  unaffected. 

5 . 5  Influence  of  Wave  Acceleration 

The  velocity  field  due  to  the  wave  motion  is,  of  course,  much  larger 
than  the  random  velocity  field  due  to  the  turbulence.  We  must  consider 
how  the  droplets  respond  to  the  wave  velocity  field,  and  how  this  affects 
their  motion. 

We  need  an  estimate  for  the  relative  velocity  between  the  droplet 
and  the  water  due  to  the  wave  acceleration.  This  may  be  obtained  from  the 
linearized  equation  for  the  droplet  motion.  We  begin  with 

dvj/dt  *  (uj  -  (5-14) 

where  v  is  the  droplet  velocity,  u  is  the  fluid  velocity  at  the  droplet 
location,  and  we  are  considering  a  horizontal  component,  the  1-component, 
aj  is  the  horizontal  droplet  time  constant;  the  horizontal  and  vertical 
time  constants  are  different,  as  a  result  of  linearization  about  a 
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terminal  velocity  which  is  in  the  non-linear  regime  of  the  drag  law.  We 
consider  a  horizontal  component  in  order  to  eliminate  the  effect  of  the 
terminal  velocity,  which  is  not  germane  to  this  discussion.  If  we  suppose 
a  monochromatic  wave  of  the  form 

n  ■  h  exp[ix(x  -  c t )  ]  (5-15) 

where  n  is  the  instantaneous  surface  elevation,  and  h  is  the  maximum  value 
of  the  surface  elevation,  then  the  potential  function  for  the  wave  motion 
is  given  by 


$  ”  i  c  h  exp[-Kz  +  ix(x  -  ct)]  (5-16) 

c  is  the  phase  velocity  of  the  wave,  and  x  its  wave  number,  z  is  measured 
positive  downward.  Substituting  the  velocity  field  resulting  from  (5-16) 
in  the  equation  (5-14),  evaluated  at  x  «  0,  we  find  that  the  particle 
velocity  vj  is  of  the  form 


vi  ■  o  exp[-xz  -  ixct) 


(5-17) 


where  o  is  given  by 


u  ■  -x  c  h/(l  -  ixcai ) 

From  this  we  may  easily  obtain  the  relative  velocity  vj  -  uj , 
mine  the  mean  square  value,  to  obtain 

((vi  -  ui )2} 1/2/{u12}1/2  .  Kcai/[i  +  (xcai)2]l/2 


(5-18) 


and  deter- 


(5-19) 
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where  we  are  using  f  }  to  indicate  an  average. 

Now,  we  must  estimate  the  values  of  the  various  quantities  appearing 
in  these  expressions.  The  horizontal  time  constant  is  given  by 

al8  “  VT  (5-20) 

where  Vj  is  the  true  (as  opposed  to  the  Stokes)  value.  This  is  given  by 
(using  equation  (5-13)) 


VTS  *  VT [ 1  ♦  (3/ 16) R]  (5-21) 

If  we  consider  our  largest  particle,  which  will  have  the  largest  time  con¬ 
stant,  we  have  VT  *  0.0816Vts»  with  R  *  60.  Using  expression  (5-6) 
for  the  Stokes  terminal  velocity,  we  have  for  Pp/pf  3  0.7,  d  *  1650 

Urn,  aj  -  3.7xlO“3  sec,  while  for  Pp/pf  3  0.95,  d  *  3000  pm,  aj  * 

—  3 

2.04x10”  sec.  We  must  now  contrast  this  with  the  wave  frequency.  The 

phase  velocity  of  the  wave  is  always  of  the  order  of  the  wind  speed  in  the 

atmosphere  Ua,  so  that  the  wave  frequency  ec  is  of  the  order  of  g/Ua.  \ 

Hence,  if  the  wind  speed  is  of  the  order  of  10  m/s,  the  wave  frequency 

will  be  of  the  order  of  1  rad/sec.  Hence,  the  quantity  that  occurs  in 

equation  (5-19) 


Kcal  "  3xlO"3  (5-22) 

Hence,  from  equation  (5-19),  the  r.ra.s.  relative  velocity  is  of  the  order 
of  0.1Z  of  the  orbital  velocity  of  the  wave  under  the  worst  conditions. 

Thus,  even  under  the  most  severe  circumstances,  the  droplet  will 
simply  follow  the  orbital  motion  of  the  wave.  Since  the  velocity  field  of 
a  wave  is  non-dispers ive  to  lowest  order,  the  dispersion  of  the  droplets 
will  be  the  result  entirely  of  the  rotational  component  of  the  velocity 
field,  the  turbulence.  We  could  formally  write  the  droplet  velocity  as 
the  sum  of  the  orbital  velocity  and  the  turbulent  velocity  (  plus  the 
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terminal  velocity  and  the  velocity  of  the  Langmuir  circulation),  but  we 
will  not  bother  to  do  this  formally,  and  will  simply  suppress  the  orbital 
velocity  in  what  follows. 

5.6  Particle  Inertia 

We  may  carry  out  a  calculation  similar  to  that  in  §  5.5  to  see 
whether  the  inertia  terras  in  the  particle  equation  are  important  for  the 
turbulent  component.  There  is  a  possibility  that  they  could  be  neglected, 
since  they  may  be  neglected  for  dust  particles  in  the  atmosphere,  for 
example  (see  Lumley,  1978).  The  relative  velocity  may  be  obtained  from 
the  relative  Reynolds  number,  expressions  (5*8)  and  (5-9).  The  ratio  of 
the  relative  velocity  to  the  turbulent  velocity  is  given  by 

i(vj  -  ux)2}1/f2/{ux2}1/2  a  (n/dlR/Rfcl/^  (5-23) 
We  have  selected  for  our  worst  case  R  =  12.5,  and  n/d  is  approximately  1 

4 

for  this  case.  Rg,  =  2x10  ,  as  a  typical  value,  and  the  1/4  power  of 
this  is  about  12.  Hence,  we  may  expect  the  left  side  of  (5-23)  to  be 
about  unity.  Thus,  the  inertia  terras  may  be  expected  to  be  quite  impor¬ 
tant  so  far  as  the  turbulence  is  concerned,  as  important  as  the  other 
terras  in  the  equation. 

Fortunately,  however,  we  are  helped  by  a  happy  circumstance.  As  is 
explained  in  Lumley  (1978),  the  direct  influence  of  the  droplet  inertia 
drops  out  of  the  expression  for  the  droplet  diffusivity  in  the  asymptotic 
limit,  so  that  the  asymptotic  diffusivity  is  dependent  only  on  the  corre¬ 
lation  of  the  fluid  velocity  seen  by  the  droplet.  The  inertia  affects 
both  the  mean  square  droplet  fluctuating  velocity,  reducing  it  below  that 
of  the  fluid  velocity  seen  by  the  droplet,  and  the  integral  time  scale  of 
the  droplet  velocity,  increasing  it  above  that  of  the  fluid  velocity  seen 
by  the  droplet  by  the  same  amount,  but  the  product  remains  the  same. 
Hence,  so  long  as  our  diffusion  tiroes  are  long  relative  to  the  integral 
t iree  scales  of  droplet  motion,  we  may  ignore  droplet  inertia  in  calcula¬ 
ting  droplet  di  f fusivities. 
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We  may  make  a  crude  estimate  of  the  droplet-attached  time  scale  by 
first  estimating  the  droplet-attached  mean  square  fluctuating  velocity. 
The  expression  for  the  mean  square  droplet  fluctuating  velocity  (in  the 
horizontal  direction)  is  given  by 

I'M2}  *  ({ui2}/ai)Jexpl-T/a;]g1,(T)dT  (5-24) 

where  the  range  of  integration  is  the  right  half-line.  We  are  again  using 
(  }  to  indicate  an  average,  and  gjj  is  the  autocorrelation  of  the  fluid 
velocity  seen  by  the  droplet.  We  know  from  section  5.2  that  the  droplet 
motion  may  be  considered  essentially  Eulerian,  so  that  gi  i  »  gCV-j-x), 
where  g(x)  is  the  Eulerian  transverse  spacial  correlation.  We  may  pick  a 
simple  form  for  this  function: 

g(x)  »  (1  -  5/4)exp{-£/2] ,  5  =  x/L22  (5-25) 

where  L22  is  the  transverse  integral  scale,  one-half  the  longitudinal 
scale  in  an  isotropic  turbulence.  Expression  (5-24)  can  be  evaluated  to 
give 

fvl2}/(vUZ}  "  Y(1  ♦  4y>/(1  ♦  2Y)2  (5-26) 

where  Y  *  L22/Vt*i  “  1*2  2  /  1 2 .  that  is,  the  ratio  of  the  transverse 
integral  scale  to  the  distance  the  droplet  moves  (at  the  terminal 
velocity)  during  one  time  scale.  The  transverse  integral  scale  is  roughly 
1/4,  or  about  0.5  meters  in  our  simple  example,  if  we  take  l  to  be  the 
wave  height,  at  2  meters.  We  have  already  estimated  a\  ■  3xl0-3  sec. 
Hence,  we  obtain  Y  *  5.67xl03.  It  is  clear  that  any  other  set  of 
reasonable  assumptions  will  not  change  the  order  of  magnitude  of  this 
quantity.  Evaluating  expression  (5-26),  we  obtain  0.999.  Hence,  the  mean 
square  fluctuating  droplet  velocity  and  the  mean  square  fluctuating  fluid 


velocity  are  essentially  the  same,  although  the  Last  ant aneous  velocities 
are  not,  due  to  the  phase  lag  induced  by  the  droplet  inertia.  Since  the 
droplet  attached  mean  square  fluctuating  velocity  is  the  same  as  the  mean 
square  fluctuating  fluid  velocity  seen  by  the  droplet,  the  droplet- 
attached  time  scale  must  be  the  same  as  the  time  scale  of  the  fluid  velo¬ 
city  seen  by  the  droplet. 

Since  the  oil  droplets  are  operating  in  the  Eulerian  regime,  we  may 
estimate  the  integral  time  scale  of  the  fluid  velocity  seen  by  the  droplet 
as  the  Eulerian  transverse  integral  scale  divided  by  the  droplet  mean  vel¬ 
ocity  relative  to  the  fluid,  the  terminal  velocity.  We  may  estimate  the 
Eulerian  transverse  integral  scale  as  about  one-quarter  of  the  wave 
height.  We  have  been  taking  as  an  example  a  ca9e  of  waves  of  two  meters 
height,  but  here  we  must  use  the  wave  heigV  -8  corresponding  to  the  various 
cases  in  tables  7-2,  7-3,  and  7-4.  Since  the  wave  slope  is  constant  at 
.088,  the  wave  heights  are  respectively  (for  wave  numbers  of  1/m.  0.5/m 
and  0.25/m)  0.088ro,  0.176m  and  0.352ra,  giving  an  Eulerian  transverse 
integral  scale  of  0.022ra,  0.044m,  and  0.088ra. 

The  time  required  to  rise  up  Che  side  of  the  Stonunel  retention  zone 
may  be  estimated  from  the  tables  in  section  7,  using  the  dimensionless 
distance  St  and  the  scaling  velocity*  Proceeding  in  this  way,  we  find 
that  the  time  required  to  rise  up  the  side  and  the  integral  time  scale  of 
the  fluid  velocity  seen  by  the  droplet  scale  the  same  way  with  the  wave 
height,  and  we  have  a  single  ratio  for  each  of  the  cases,  giving  a  ratio 
of  the  integral  time  scale  to  the  convection  time  of  0.15,  0.2  and  0.26  in 
the  cases  A,  B  and  C.  Hence,  we  see  that  in  the  worst  case  the  rising 
droplet  spends  about  four  integral  time  scales  rising  up  the  side,  so  that 
the  diffusion  coefficient  is  in  the  asymptotic  regime  long  before  the 
droplet  reaches  the  top.  Hence,  it  is  legitimate  to  use  the  asymptotic 
form  of  the  diffusivity  for  our  problem,  and  we  may  ignore  droplet 
inert ' a . 
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5 . 7  The  Crossing-Trajectories  effect  and  the  Pi  f fusivit ies 

The  construction  of  the  dif fusivities  within  our  restrictions  is 
described  in  Lottey  (1979).  Briefly,  the  argument  is  this:  Consider  the 
vertical  velocity.  For  large  Reynolds  number  flow9,  such  as  the  wave- 
associated  turbulence  under  consideration,  the  large  scale  structure  of 
the  velocity  field  may  be  considered  to  be  determined  by  the  integral 
length  scale  and  the  r.ra.s.  velocity.  The  droplet  velocity  autocorrela¬ 
tion  will  be  a  function  of  the  Eulerian  longitudinal  integral  length 
scale,  the  r.ra.s.  velocity,  the  droplet  terminal  velocity  and  the  time 
lag.  These  form  two  dimensionless  groupings,  one  involving  the  r.ra.s. 
fluctuating  velocity,  and  the  other  the  terminal  velocity.  If  the  termin¬ 
al  velocity  is  small  compared  to  the  r.ra.s.  velocity,  the  case  is  essen¬ 
tially  Lagrangian,  while  if  the  reverse  is  true,  the  case  is  essentially 
Eulerian.  If  the  correlations  in  these  two  limiting  cases  are  assumed  to 
have  the  same  shape,  then  the  function  of  two  variables  can  be  replaced  by 
a  function  of  a  single  variable  -  essentially  the  isocorrelation  contours 
are  assumed  to  have  the  form  of  ellipses.  There  is  a  constant  relating 
the  two  axes  which  is  determined  by  an  integral  constraint,  assuring  that 
the  integrals  along  the  two  axes  give  the  correspond  ing  (Eulerian  and 
Lagrangian)  integral  scales.  Finally,  the  integral  may  be  carried  out  to 
give  the  asymptotic  dispersion,  giving 

I>2  ■  ( u£ / 3 )  1 1  ♦  (4/9) (V-j/u) 2  (5—27) 

independent  of  the  form  of  the  correlation.  Similar  arguments  may  be  made 
for  the  horizontal  direction,  but  here  a  difficulty  arises.  The  Eulerian 
transverse  correlation  goes  negative  due  to  the  requirement  that  the  total 
mass  flux  across  any  plane  must  be  zero.  However,  the  Lagrangian  correla¬ 
tion  does  not  go  negative;  hence,  it  is  not  possible  to  assume  that  the 
form  of  the  correlation  in  the  two  limiting  situations  is  the  same.  How¬ 
ever,  measurements  by  Snyder  4  Luraley  (1971)  indicate  that  the  droplet- 
attached  velocity  autocorrelations  and  the  Eulerian  transverse  correlation 
are  similar  down  to  very  small  values  of  the  terminal  velocity.  Hence, 
proceeding  in  the  same  way,  we  obtain 

D1  "  d2^2  (5-28) 
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This  expression  for  Dj  does  not  give  the  right  value  as  V-p  ♦  0,  since  it 
gives  half  the  Lagrangian  diffusivity.  This  is  due  to  the  fact  that  the 
droplet-attached  correlation  changes  rapidly  as  the  Lagrangian  case  is 
approached,  and  no  adequate  approximat ion  is  known.  However,  down  to  very 
small  values  of  Vj/u  it  is  expected  that  the  expression  given  for  Dj 
will  be  adequate. 

5 .8  The  Fluid  Diffusivity 

The  quantity  appearing  in  (5-27),  uA/3,  is  the  fluid  diffusivity. 
We  need  an  estimate  for  this.  To  obtain  such  an  estimate,  we  have  anal¬ 
yzed  measurements  taken  in  Lake  Ontario  by  Donelan,  and  reported  in  Kitai- 
gorodskii  et  al  (1981).  These  measurements  were  taken  under  waves  from  5 
to  15  cm  in  height,  having  wavelengths  between  three  and  ten  meters,  with 
winds  (at  ten  meters  height  above  the  mean  water  surface)  between  five  and 
fifteen  meters  per  second.  These  conditions  are  quite  typical  of  situa¬ 
tions  of  interest  here  with  the  exception  of  the  wave  height,  which  we 
would  expect  to  be  larger  in  the  open  ocean,  with  either  longer  fetches, 
or  greater  development  time.  The  measurements  have  been  parameterized  as 

u/u^  *  (0 .28(hm/u*)  ♦  2. 30]exp[-<z/2)  (5-29) 

where  the  numerical  coefficients  are  independent  of  sea  state.  The 
wavenumber  of  the  exponential  decrease  is  not  as  certain,  the  factor  rang¬ 
ing  from  about  0.4  to  0.9;  we  have  selected  0.5  as  a  compromise  value 
because  it  results  in  a  considerable  simplification,  as  we  shall  see 
shortly . 

Unfortunately,  the  measurements  of  the  dissipation  rate  £  of  turbu¬ 
lent  kinetic  energy  are  much  less  reliable,  and  little  can  be  learned  from 
them,  other  than  that  the  length  scales  predicted  are  of  the  order  of  the 
wave  height.  In  this  connection,  we  may  mention  the  estimates  of  £  sug¬ 
gested  by  Milgram  et  al  (1978),  who  propose  10g2/ui  for  breaking  waves.  If 
we  calculate  from  this  the  logarithmic  decrement  of  the  wave  in  periods, 
we  find  that  the  wave  will  decrease  in  amplitude  by  1/e  in  1/10  of  a 
period.  This  is  clearly  much  too  large,  and  the  values  of  n  they  calcu¬ 
late  using  this  value  of  £  are  much  too  small.  Clearly  breaking  waves 
will  have  larger  dissipation  than  non-breaking  ones,  but  one  would  expect 
the  loss  of  energy  to  the  breaking  to  take  a  number  of  periods. 
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Rather  than  use  unreliable  experimental  results,  we  have  used  a 
theoretical  model  to  obtain  the  values  of  e,  or  equivalently,  a  length 
scale.  If  we  suppose  that  the  Reynolds  stress  induced  in  the  turbulence 
in  the  wave  by  the  strain  rate  of  the  wave  motion  is  related  to  the  latter 
by  an  eddy  viscosity,  as 

-{ujujl  +  (q2/3)«ij  -  vTS{ j  (5-30) 

where  S^j  is  the  strain  rate  of  the  oscillatory  motion  of  the  wave,  then 
the  rate  of  production  of  turbulent  energy  is  given  by  vT^ij^ij» 
which  can  be  written  as  u^/4  (Tennekes  &  Luraley,  1971).  Writing  \>t  as 
uA/3,  we  obtain  (computing  the  oscillatory  motion  of  the  wave) 

*u/3  -  vT  *  u2//6hc<2exp[-»cz]  (5-31) 

Using  the  behavior  we  assumed  for  u,  we  obtain  in  this  way  a  v-j-  which  is 
independent  of  depth  of 

vT  =  u*2[0.28(hu>/uA)  +  2. 30]  2//6ho»ic  (5-32) 

and  we  have  used  this  expression  in  the  form  for  Dj  and  D2. 
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A  COMPUTATIONAL  METHOD  TO  IMPLEMENT  THE  MODEL 


Before  explaining  procedures  used  to  solve  the  diffusion  equation 
(4-11),  we  review  the  features  that  have  already  been  discovered,  and  give 
a  descriptive  account  of  the  physical  processes  taking  place. 

First,  we  have  noted  that  diffusion,  although  an  essential  process, 
is  weak.  If  ignored,  particles  move  on  lines  of  constant  T  as  indicated 
in  Figure  4-4  .  If  an  oil  particle  is  placed  on  a  ¥  line  outside  a 
retention  zone,  then  it  will  rise  along  that  line  until  it  reaches  the 
surface,  z  =  0.  Since  all  particles  reaching  the  surface  reside  in  a 
sheet  of  zero  thickness,  the  subsequent  motion  of  one  of  these  particles 
cannot  be  inferred  from  the  figures.  Nevertheless,  the  behaviour  of  a 
particle  on  the  surface  is  clear.  According  to  the  diffusion-free 
approximation,  it  must  stay  on  the  surface  and  move  towards  LC  convergence 
lines.  Thus,  in  this  simplified  view,  the  entire  volume  outside  of  the 
retention  zones  is  swept  free  of  oil,  all  of  the  oil  initially  there  being 
transported  to  the  surface  and  then  swept  into  windrows.  The  final 
picture  would  be  one  where  all  of  the  oil  is  concentrated  in  a  line  on  the 
surface  -  the  windrow  -  and  possibly  oil  trapped  below  the  surface  in  a 
SRZ . 

With  an  infinitesimal  amount  of  diffusion  acting,  any  oil  trapped  in 
a  SRZ  would  ultimately  be  uniformly  distributed  throughout  the  SRZ.  This 
is  one  effect  of  the  diffusivity.  Since  the  simplified  view  gives  rise  to 
a  singular  distribution  of  surface  oil,  with  infinite  concentration  grad¬ 
ients,  we  know  that  a  second  effect  of  small  but  non-zero  diffusivities  is 
to  spread  the  surface  oil  out  into  a  band  of  non-zero  thickness  and  non¬ 
zero  width  in  the  windrow.  A  third  effect  is  the  slow  transport  of  oil 
across  the  boundary  of  the  SRZ.  Since  the  diffusivities  are  weak,  oil 
droplets  transferred  across  the  side  of  an  SRZ  form  a  thin  layer  adjacent 
to  the  retention  zone  boundary;  once  outside  they  rise  up  towards  the  sur¬ 
face  within  the  thin  boundary  layer.  This  is  depicted  in  the  sketch  of 
Figure  6-1. 
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Figure  6-1.  -Sketch  of  SRZ  ;  its  windrow  supply  source,  and 

boundary  layer  leakage  rising  alon^  the  sides. 

If  there  were  no  other  effects  of  diffusion,  the  boundary  leakage 
would  eventually  empty  the  retention  zone,  and  alt  oil  would  be  collected 
into  surface  windrows  with  vertical  and  horizontal  spreads  that  are 
functions  of  V^,  the  Langmuir  circulation  velocity  scale  V,  and  the 
dif fusivities.  If,  however,  the  top  of  the  SRZ,  defined  as  the  first 
point  below  the  windrow  at  which  w  *  -  V'j,  lies  within  a  distance  below 
the  surface  characteristic  of  the  vertical  extent  of  the  windrow,  then  the 
vertical  concentration  gradient  in  the  windrow  can  produce  a  vertical  flux 
of  oil  particles  back  into  the  retention  zone.  Under  these  circumstances 
it  is  possible  to  arrive  at  a  steady  balance  between  the  loss  of  oil 
across  the  sides  of  the  SRZ  and  the  injection  of  oil  from  the  windrow  into 
the  SRZ. 

The  purpose  of  this  section  is  to  compute  the  equilibrium  balance 
that  ultimately  develops.  Since  the  dif fusivit ies  are  so  small,  a 
straightforward  calculation  (e.g.,by  a  finite  difference  algorithm)  of  the 
field  as  a  whole  is  not  feasible.  A  mesh  fine  enough  to  resolve  the 
various  thin  regions  would,  if  applied  to  the  field  as  a  whole,  result  m 
excessive  storage  requirements  and  computing  time.  Consequently,  the 
calculation  is  diviaed  into  parts  that  take  advantage  of  the  problem 
structure . 

The  boundary  layer  on  the  rising  side  of  the  zone  is  calculated 
analytically;  the  onLy  parameter  that  enters  the  finnl  result  is  an 


effective  diffusion  time  along  the  rising  side  and  involving  the 
horizontal  and  vertical  d i f fusivi t ies  (reflecting  the  change  of  direction 
of  the  droplet  motion  as  it  rises)  as  well  as  the  mean  droplet  speed.  The 
analytical  calculation  can  be  done  once  for  all  cases.  To  apply  the 
result,  the  effective  diffusion  time  must  be  determined  for  each  Langmuir 
cell  and  each  set  of  droplet  characteristics. 

The  region  between  the  top  of  the  Stommel  retention  zone  and  the 
surface  must  be  solved  numerically.  The  dispersion  in  this  region  is 
dependent  upon  two  parameters:  the  dimensionless  distance  of  the  top  of 
the  zone  from  the  surface,  and  the  thickness  of  the  droplet  boundary  layer 
entering  from  the  side.  The  calculation  is  delicate  because,  even  in  this 
reduced  region,  the  gradients  are  large,  and  it  is  necessary  to  write  the 
program  in  such  a  way  that  the  droplets  are  conserved.  From  such 
calculations  we  determine  the  relative  concentration  within  the  Stomrael 
retention  zone  and  in  the  windrow  for  given  sea  state,  Langmuir  cell  and 
droplet  parameters. 

Under  other  support,  measurements  of  turbulent  intensities  and 
scales  under  various  sea  states  have  been  analyzed,  and  computational 
modeling  of  the  turbulence  in  the  surface  mixed  layer  of  the  ocean  in  the 
presence  of  waves  has  been  carried  out.  By  a  combination  of  these  two,  we 
have  arrived  at  a  distribution  of  turbulent  droplet  di f fusivit ies  as  a 
function  of  sea  state  as  described  in  the  previous  section.  The.  o 
diffusivit ies  enter  into  the  determination  of  the  effective  dimensionless 
distance  of  the  top  of  the  Stommel  retention  zone  from  the  surface,  and  of 
the  effective  diffusion  time  for  the  rising  droplets  on  the  edge  of  the 
Zone. 

6. 1  Preliminary  Survey 

We  now  put  the  picture  outlined  above  int  >  mathematical  form. 
Consider  the  diffusion  equation  in  its  dimensionless  form  (4-17).  It  is 
convenient  to  choose  D2(0)  to  serve  as  a  scale  for  the  di ffusivities,  and 
to  write 

Aj (z)  -  D!(z)/D2(0),  A2(z)  "  D2(z)/D2(0),  e  -  k0D2(0)/V; 

( 6*1 ) 

thus  A  2(0)  *  1  and  £  is  a  small  parameter.  In  this  notation,  the 
diffusion  equation  and  boundary  conditions  on  C(y,z)  for  steady  conditions 
are  expressed  as 
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vC  +  (w  +  V'T)C,Z  =  e{(iiC  ),y  ♦  (A2C,z),z} 

J  (a) 

C,y(y,0)  -  V'x  C(y ,0)  (b) 

C,y(0,z)  *  C,y(Y2 ,z)  *  0  (c) 

C  +  0  as  2  ♦  -  “  (j) 

where  all  quantities  are  dimensionless,  V’j  *  V^/V;  y  =  0  is  assumed 
to  be  the  location  of  a  downwelling  plane,  and  y  *  Y2  the  dimensionless 


location  of  an  upwelling  plane,  so  a  single  cell  in  the  given  Langmuir 
circulation  system  lies  in  0  y  <.  Y2  .  Furthermore,  an  SRZ  lies  below  the 
LC  windrow.  Assuming  symmetry  in  y,  half  of  the  SRZ  lies  in  the  LC  cell 
considered,  and  its  boundary  f  =  0  includes  an  interval  on  the  z-axis  from 
z  *  z_,  the  lower  boundary  of  the  SRZ  to  z  *  z+ ,  the  top  of  the  SRZ. 
The  sketch  of  Figure  6-2  shows  the  situation  analyzed. 


Figure  6-2.  Definition  sketch  for  diffusion  calculation. 
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Consider  the  computing  requirements  for  a  finite  difference  analysis 
of  this  problem  for  e  *  0.005,  a  typical  value  for  conditions  to  be  dealt 
with;  and  Y2  *  2,  which,  upon  referring  to  the  Figures  in  $4,  is  one  of 
the  cases  (1)  we  analyse.  If  we  suppose  that  we  were  to  compute  with  a 
raesh  of  Ay  »  A z  m  0.1  and  imagine  iterating  or,  equivalently,  calculating 
an  initial  value  problem  until  a  convergence  to  a  steady  state  is  reached, 
then  this  mesh  is  about  the  largest  we  can  use  to  avoid  unacceptable 
levels  of  artificial  diffusion.  With  this  mesh,  the  computational  domain 
0<y<2,  0<-z<5  (this  is  the  maximum  penetration  of  the  LC)  would  be 
resolved  with  a  20x50  grid, or  1000  mesh  points.  To  achieve  a  steady 
state,  we  must  integrate  an  initial  value  problem  for  a  time  comparable  to 
the  diffusion  time  in  the  problem.  Based  upon  the  horizontal  distance  Y2 
(which  may  be  optimistic),  since  the  vertical  extent  of  the  first  is 
larger  than  Y2 ,  this  is  a  dimensionless  time  t  *  Y2 2/e  *  ,8xl03.  Numeri¬ 
cal  stability  requirements  demand  that  the  time  step  be  restricted  to  the 
smaller  of  the  advective  time  scale  based  upon  the  spatial  mesh  or  the 
diffusion  time  based  upon  the  mesh.  The  former  is  operative  here,  and  the 
time  step  is  on  the  order  of  0.1,  Thus,  a  minimal  calculation  at  field 
sight  requires  approximately  10,000  computations  at  each  1000  mesh  points, 
leading  to  more  than  108  arithmetical  operations.  While  this  is  manage¬ 
able,  large  concentration  gradients  in  the  windrow  region  require  a  sub¬ 
stantially  finer  grid  in  that  critical  area,  and  this  we  estimate  would 
double  the  storage  required  and  increase  the  overall  number  of  operations 
by  a  factor  exceeding  105. 

While  better  numerical  methods  are  probably  possible  that  can  reduce 
the  computational  effort  considerably,  it  seemed  to  us  that  the  end  result 
would  lead  inevitably  to  very  large  computational  problem.  Rather  than 
fight  the  smallness  of  e  we  preferred  to  choose  a  method  that  exploits 
it.  This  we  have  done  by  dividing  the  Langmuir  ceil  into  the  regions  de¬ 
scribed  in  the  introduction  to  this  section.  These  regions  are  schemati¬ 
cally  indicated  in  Figure  6-3;  the  decomposition  of  the  problem  and  the 
thickness  of  the  various  regions,  are  the  subject  of  this  subsection. 
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Figure  6-3.  Decomposition  of  the  domain. 

We  begin  by  noting  that  in  (6-1)  v  and  w  are  0(e°).  Regions  I  and 
II  are  those  defined  by  the  limit  problem,  e  *  0,  and  consists  of  the  SRZ 
(I)  and  the  region  outside  it  (II).  Since  the  SRZ  has  a  uniform  oil 
concentration,  and  since  the  solution  of  the  steady  problem  is,  as  we  have 
already  discussed,  invariant  with  respect  to  an  arbitrary  multiplication 
factor,  we  take  c  s  1  in  the  SRZ.  In  region  II,  being  free  of  oil  if  €  ■ 
0,  c  *  0.  Regions  I  and  II  are  the  "outer”  regions  (van  Dyke,  1975)  and 
are  separated  by  an  "inner"  transition  region  III  in  which  oil  from  inside 
the  SRZ  escapes  and  travels  up  towards  the  surface.  This  transition 
region  has  a  thickness  proportional  to  /e,  so  diffusion  tangent  to  the  SRZ 
boundary  is  negligible  compared  to  that  normal  to  it,  and  is  therefore 
ignored.  This  approximation  cannot  be  made  to  satisfy  the  boundary 
condition  at  y  *  0;  thus  region  IV,  in  which  the  full  equations  mist  be 
solved,  is  needed.  Region  IV,  as  will  be  shown,  is  a  layer  along  the  y 
axi9  with  thickness  of  order  le | z+| /V'xl ^ ^ . 

6.2  Region  III.  Escape  of  Oil  from  the  Retention  Zone. 

To  treat  the  transfer  of  oil  across  the  curved  boundary  of  the  SRZ, 
we  adopt  intrinsic  coordinates  lying  along  and  normal  to  the  droplet 
pathlines.  The  droplet  velocity  and  "st reamfunct  ion"  are  (u,  w  +  Vf) 
and  V ,  respectively,  and  Y  *  0  is  taken  to  be  the  bou  ‘ary  of  the  SRZ. 
Let  s  be  arc  length  measured  along  the  SRZ  boundary  beginning  at  the  bottom 
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(z  “  z_)  and  let  n  be  Che  coordinate  normal  to  the  boundary  and  directed 
towards  the  interior  of  the  retention  zone  with  n  “  0  coinciding  with  V  ■ 
0.  Let 

q  “  [v2  ♦  (w  +  V'j)2)1/2  (6-2) 

be  the  speed  of  particles  as  they  rise:  this  depends  upon  s  and  n,  but  the 
difference  of  q  (and  4j  or  &2  aa  well)  from  its  value  at  n  *  0  to  any 
other  point  in  region  III  is  higher  order  (0  (/e)  and  can  be  neglected. 
With  the  choices  made 


I1  *  /gq(s  »n^n  *  <5  n 


(6-3) 


to  the  order  considered. 

If  we  introduce  a  boundary  layer  stretch  in  (6-1)  by  setting 

n  -  /e  N,  f  -  /e  n  (6-4) 

(where  n  *  q  N),  equations  (6-1)  take  the  form 

-  [4i(n,y)2  +  4z(n,l)2]C,NN  (6-5) 

accurate  to  0(/e).  This  equation  is  more  conveniently  expressed  in  terms 
of  the  stretched  streamfunction  n,  using  (6-4);  it  then  takes  the  form 

qC,s  -  [4;  (w  +  V'x)2  +  4zv2]C>T,n.  (6-6) 

Finally,  we  note  that  &i,  4Z ,  v,  w,  q  all  may  be  considered  to  vary  only 
with  s  in  III,  and  so  we  make  a  transformation  to  a  new  coordinate 


K 


/ q“ 1 [ Ai (w 


o 


V'T)Z  +  42V2 Ids 


(6-7) 


to  arrive  at  the  form  of  the  concentration  equation  that  we  deal  with  in 
region  III: 
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“  c.nn- 


(6-8) 


This  is  the  heat  equation  in  its  simplest  form. 

The  boundary  conditions  on  (6-8)  are  matching  with  regions  X  and  II, 
or 


C  *  1  S!  II  »  '  * 
C  *  0  as  n  » 


(6-9) 


and  it  is  assumed  that  region  HI  starts  at  (  =  0,  the  bottom  of  the  SRZ, 
at  which  point  C  3  1  for  n  >  0  and  C  =  0  for  n  <  0. 

The  solution  for  region  III  is 

c  -  (l/2)erfc(-n/2/0  (6-10) 


where  erfc  is  the  complementary  error  function  (Abramowitz  and  Stegun, 
1964).  As  indicated  earlier  the  solution  (6-10)  cannot  be  made  to  satisfy 
the  boundary  condition  (6-1  c),  and  we  are  therefore  forced  to  consider 
region  IV. 

6 .3  Region  IV.  Oil  Transfer  from  Windrow  to  Retention  Zone. 

In  order  to  meet  the  conditions  imposed  at  y  *  0,  the  end  of  the 
rising  droplet  layer,  diffusion  in  the  direction  along  the  retention  zone 
boundary  must  be  restored.  Since  this  failure  of  the  boundary  layer 
solution  occurs  at  the  top  of  the  SR2,  at  z+  and  y  *  0,  directions  of  s 
increasing  and  n  increasing  are  anti-parallel  to  y  increasing  and  z 
increasing.  Clearly,  the  additional  term  is  important  because  gradients 
in  the  y-direction  become  as  large  as  those  in  the  z  (i.e.,  n)  direction 
in  the  neighborhood  of  y  *  0.  To  make  this  explicit,  a  stretching  in  the 
y-direction  with  the  same  magnif ication  provided  in  the  n-direction  is 
required. 

To  begin  the  analysis  of  this  region,  suppose  we  try  a  stretch 
e"a,  a  >  0,  so 

y  *  Y’Ea  (6-n) 
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We  continue  to  look  at  the  boundary  layer,  so  we  retain  the 
stretch  in  the  n-direction,  but  re-written  in  terms  of  2.  In  particular, 
let  us  put 

2  -  2+  +  Z'A  (6-12) 

so  that  Y  *  0,  Z  ■  0  is  the  top  of  the  retention  eone.  Substitute  (6-11) 
and  (6-12)  into  (6-1);  the  dominant  terms  are 

e-a  v(ea  y',z+  +  Jez )C>y,  +  l„(ea  V  , 

z+  ♦  /e  Z)+V'T]C,Z</A  (6-13) 

*  e  (j,+  +  /e  Z'  )C,y'y*  *  *2^z+  +  2’ )C,z'Z' . 

Now  the  velocity  field  may  be  expanded  in  a  series.  Since  v(0,  z)  Z  0, 
v(e»  V  ,  *♦  ♦  A  Z’)  -  e8  Y'v,y(0,z*)  ♦  0(e1/2+8). 

At  the  top  of  the  retention  zone, 

w(0,  z+)  ♦  V'j  “  0 

hence  the  function 

wU8Y',z+  ♦  Az')  ♦  V'T  -  AZ'w,x(0,z+) 

since 

w,y  (0,  z+)  •  0,  and  from  continuity 
w,z(0,z+)  “  -  V,y(0,Z+) 

and  we  let  a  be  the  common  value.  In  terms  of  the  strain  rate  a,  then, 
(6-13)  is 

-u  Y'C.y'  ♦  a  Z'C,Z'  »  e1_286i  (z+)C,y'y'  *  d2(z+>C,z'Z'  (6-14) 
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and  we  choose  a 
restore. 


1/2  to  retain  the  diffusion  term  we  had  intended  to 


Equation  (6-14)  can  be  written  in  a  neater  form  by  taking 

Y'  -  Y(A1(z+)/a)l/2,  Z'  -  Z(A2 (z+)/a)l/2  (6-15) 

in  which  case  (6-14)  is 

-Y  C,y  +  Z  C,z  -  C,yy  +  C  ,22  (6-16) 

and  this  is  the  equation  that  we  solve  in  region  IV. 

Retracing  the  transformations  used,  we  see  that  final  coordinate 
choice  is  related  to  the  y  and  z  coordinates  as  follows 

y  “  YUAj/a)!/2,  z  =  z+  +  ZUAj/a)1/2  (6-17) 

Note  that  (cAj/a)!/2  is  the  (dimensionless)  thickness  of  a  boundary 
layer  formed  by  the  balance  between  (turblent)  diffusion  and  mean  strain 
rate. 

Before  we  proceed  to  solve  the  problem  in  region  IV,  we  make  the 
following  important  observation.  If  region  IV  does  not  overlap  with  the 
windrow,  then  all  oil  entering  it  from  the  rising  boundary  layer  will 
leave  through  the  top  of  region  IV  and  rise  to  the  surface;  thus  the  SRZ 
will  be  swept  free  of  oil  unless  there  is  an  overlap.  Thus,  we  mjst 
require  that  the  top  of  the  SRZ  have  (dimensionless)  depth 

d  •  jz+|  *  (Kea/Az)1/2 

if  there  is  to  be  oil  in  the  SRZ.  Since  |z+|  -  0(/e)  we  may  replace 
Aj(z+)  by  Aj  (0) ,  A2(z+)  *  A2(0)  “  1,  and  we  may  also  relate  the  strain 
rate  o  to  V'y.  In  particular,  since 

w  +  V'T  -  o  (z  -  z+)  -  o  (z  +  d) ,  (6-18) 

if  z  “  0  is  included  in  IV, 

v't  “  o  d,  (6-19) 


since  w  ■  0  at  z  *  0. 

Assuming  that  the  overlap  exists,  we  let 
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z  -  Y 


(6-20) 


locate  the  free  surface  z  “  0,  so 

Y  -  d(a/eA2)!/2  *  dCa/e)1/2  (6-21) 

on  recalling  that  d2(0)  «  1. 

The  boundary  condition  at  z  “  0  now  takes  the  form 

ClZ«  Y  C  at  z  «  Y  (6-22) 

The  problem  in  region  IV  may  now  be  summarized.  We  solve  (6-16) 
subject  to  the  boundary  condition  (6-22);  the  condition 

C »Y  ■  0  at  Y  -  0  (6-23) 

arising  from  (6-1);  matching  with  the  SR2,  or 

C  ♦  1  as  Z  ♦  -  »  ;  (6-24) 

and  matching  with  the  rising  boundary  layer.  This  latter  condition  is  now 
derived . 

The  matching  condition  is  derived  (van  Dyke,  1975)  by  rewriting  the 
"outer"  (rising  boundary  layer)  solution  in  inner  (Y,  Z)  variables, 
expanding  for  small  e  and  regrouping.  The  outer  solution  is  given  by 
(6-10).  Since 

t  -  +  0(/e), 

where  is  the  value  obtained  by  placing  s  (6-7)  equal  to  the  total 

arclength  around  one  side  of  the  SRZ,  and  since,  from  (6-4) 

n  *  f//e 

substitution  into  (6-10)  yields 
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C  -  (l/2)erfc  [-  t/(eet)l/2j .  (6.25)  I 

We  know  that  | 

V(y  ,z)  *  -  /y(w  +  V-p)dy 
0 

Substituting  from  (6-17)  for  y  and  z  and  using  (6-18)  we  have 

C  -  1/2  erfcl-YZleAj/^t)1''2]  .  (6-26) 

Which  ia  to  hold  for  Y>1 . 

We  solve  the  problem  in  region  IV  numerically  by  a  finite  difference 
algorithm,  described  below.  In  implementing  (6-26),  at  the  right-hand 
edge  of  IV,  we  set  Y  ■  50;  this  is  large  enough  to  ensure  that  the  func¬ 
tion  (6-26) is  essentially  zero  at  the  top  boundary  for  all  interesting 
values  of  Y  and  £t. 

Note  that  the  dimensionless  quantity  eAj /£t  depends  only  on  the 
value  of  the  terminal  velocity  and  the  Langmuir  cell  velocity  and  length 
scales.  This  variable  determines  the  width  of  the  incoming  concentration  \ 

distribution  at  the  right  edge  of  the  field  under  the  windrow.  We  comput¬ 
ed  the  values  of  this  variable  for  a  number  of  different  configurations. 

6.4  Computational  Procedure 

We  computed  solutions  of  equation  (6-5)  in  a  region  extending  from 
Z  *  -2  to  Z  *  Y  >  and  from  Y  ■  0  to  Y  ■  50,  for  a  number  of  different 
values  of  Y  and  a  number  of  different  Langmuir  circulation  fields  and  ter¬ 
minal  velocities.  Equation  (6-5)  was  written  as  an  unsteady  equation 
(simply  including  the  time  derivative),  and  a  simple  forward-time  center- 
ed-space  differencing  scheme  was  used.  The  grid  points  were  centered  in 
the  cells,  and  the  equation  and  boundary  conditions  were  written  in  a  con¬ 
servative  form.  Each  computation  resulted  in  a  value  of  the  ratio  of  the 
concentration  in  the  center  of  the  windrow  (Y  *  0,  Z  *  y)  to  that  in  the 
Storamel  retention  zone.  A  typical  iso-concentration  plot  is  shown  in  Fig¬ 
ure  6-4.  The  iso-concentration  contours  become  too  close  together  in  the 
center  of  the  windrow  for  the  plotting  routine  to  display  them. 
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Figure  6-4.  Example  of  a  contour  of  oil  concentration  in  region  IV.  Upper  edge  is  the  water  surface; 

top  of  the  Stommel  retention  zone  Is  half-way  up.  Contours  become  too  dense  for  the  plot¬ 
ter  to  f  illow  in  the  upper  left  corner  (center  of  windrow).  Contour  intervals  are  0.125. 
Vertical  scale  is  expanded  by  a  factor  of  five. 


RESULTS  OF  DISPERSION  COMPUTATIONS 


The  results  of  the  concentration  computations  that  we  have  made 
using  the  method  of  §6  are  summarized  in  Tables  7-1  to  7-4. 

One  of  the  principal  measures  of  oil  below  the  surface  in  this 
analysis  is  the  ratio  of  the  maximum  surface  concentration,  which  appears 
in  the  center  of  a  windrow,  to  the  concentration  of  oil  distributed 
throughout  the  Storamei  retention  zone.  This  ratio  we  designate  to  be  C*. 
As  the  dispersion  increases,  C*  decreases.  Tables  7-1  through  7-4  show 
the  effects  of  wind  speed,  sea  state,  and  depth  of  penetration  upon  C*. 
In  these  Tables,  the  wind  speed  is  identified  by  the  first  letter  of  the 
code  (A+5  ra  s"1  ,  B+10  m  s"1  ,  C+15  ra  s-1),  the  sea  state  by  the  first  num¬ 
eral  ( 1  '*•<0  m  1  m“l,  *  *5  m“l,  3+kq  .25  m”1),  aspect  ratio  of  the  con¬ 
vective  cells  (penetration  depth/windrow  separation)  is  indicated  by  the 
last  pair  of  symbols  which  refer  to  the  dimensionless  development  time  of 
the  cells  (Tl+35  units,  T2+50  units,  T3+95  units). 

Table  7-1  displays  these  results  for  a  wind  speed  of  5  m/s  for  par¬ 
ticle  terminal  velocities  of  .26  cm/s,  .65  cm/s  and  1.3  cm/s.  The  value 
ol  C*  for  a  given  wind  speed  depends  strongly  upon  V^,  but  only  weakly 
upon  the  sea  state  or  penetration  depth  of  the  cells.  Cases  AlTl  and  A3T1 
represent  situations  occuring  for  the  same  wind  speed  and  cell  aspect  ra¬ 
tio  but  sea  states  with  wave  heights  differing  by  a  factor  of  four.  As 
can  be  seen,  the  values  of  C*  are  not  much  different.  The  last  three  rows 
compare  C*  for  the  same  terminal  velocity,  but  different  aspect  ratios 
and/or  sea  states.  Again,  the  values  of  C*  are  comparable. 

We  may  conclude  from  Table  7-1  that  C*  is  a  rapidly  increasing  func¬ 
tion  of  V^/V,  where  V  is  the  LC  convective  velocity  scale,  and  depends 
only  weakly  upon  the  sea  state  or  cell  aspect  ratio.  The  fact  that  C*  is 
not  a  large  number  means  that  for  the  cases  computed,  at  least,  vertical 
dispersion  is  effective,  and  a  considerable  amount  of  oil  is  held  in  sus¬ 
pension. 
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Table  7-1  -Concentration  Ratio  Case  A.  Wind  speed  5  m/s.  Case 
codes  as  described  in  text.  C*  is  the  ratio  of  the  maximum 
concentration  in  the  windrow  to  the  concentration  in  the  Stommel 
retention  zone.  More  extensive  data  is  reported  for  Vj  =  .0065  to 
show  that  C*  depends  only  weakly  on  <0  (represented  by  first  digit 
of  case  code)  or  upon  cell  aspect  ratio. 


Case 

vT 

m/s 

— 

c* 

A1T1 

.0026 

1.6 

.0065 

3.3 

.013 

11.0 

A3T1 

.0026 

1.6 

.0065 

3.8 

.013 

11.7 

A1T2 

.0065 

3.6 

A1T3 

.0065 

3.1 

A2T1 

.0065 

3.6 

The  main  results  of  the  dispersion  calculations  are  shown  in  Tables 
7-2  through  7-4,  for  terminal  velocities  of  0.5,  1.0,  and  0.2  cm/s, 
and  for  representative  values  of  cell  aspect  ratios  and/or  sea  states. 
These  Tables  contain  all  of  the  relevant  data  used  in  the  calculations 
(surface  value  of  di f fus ivi t ies ,  Dj(0);  e;  y)  as  well  as  C*.  In 
addition,  the  final  two  columns  of  each  Table  give  the  fraction  of  oil 
held  in  suspension  in  the  retention  zones,  and  the  depth  of  the  deepest 
retention  zone.  The  oil  fraction  is  the  percentage  of  the  total  amount  of 
oil  present  that  is  located  in  the  retention  zone,  excluding  oil  in  the 
region  IV  overlap  with  the  windrow.  These  two  entries  are  the  key  results 
of  the  dispersion  model. 

Tables  7-3  and  7-4  are  abbreviated,  and  show  only  the 
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Table  7-3  -Vy  ■  1x10"  m/ s .  Results  of  dispersion  analysis. 


Case 

Tl 

d2(0) 

m  /  s 

e 

Y 

| 

C*  1 

suspend¬ 

ed 

oil,  l 

.  _  i 

D«pth 

of 

SRZ  m 

A1T1 

Va=  3m/ a 
1/kq-  1  m 

0.00297 

0.030 

0.17 

0.619 

4.53  1 

9 

1.6 

A3T1 

Va=5  ra/s 
l/<0 *  4  m 

0.00987 

0.025 

0.17 

0.679 

5.21 

16 

6.2 

cm 

Va=15m/a 

l/Kg-  1  m 

0.00557 

0.019 

0.19 

0.257 

1.86 

24 

2.1 

C3T1 

V  a“ 1 5m/ a 
l/<0  4  m 

0.01428 

1 

0.012 

0.19 

0.321 

1 

2.19 

66 

8.4 

Tabl< 

Case 

Tl 

!  7-4  -V-] 

D2(0) 

m  /« 

r  -  2x10"' 

€ 

!  m/s .  R< 

5t 

»8Ult8  Of 

V 

dispersic 

C* 

>n  analysi 

suspend¬ 

ed 

oil,  X 

LS  . 

Depth 

of 

SRZ  m 

AlTl 

Va*5  m/s 

1/ko*  1  TO 

0.00297 

0.030 

0.17 

0.122 

1.28 

21 

2.2 

A3T1 

Va»5  m/s 
1/cq  4  m 

0.00987 

0.025 

0.17 

0.134 

1.32 

93 

9.0 

cases  corresponding  to  the  smallest  cell  aspect  ratio  computed,  since  the 
trends  are  clear  from  Table  7-2.  Table  7-4  also  includes  only  the  lowest 
sea  state;  for  the  small  value  of  Vj  used  for  this  table,  our  numerical 
results  were  felt  to  be  too  uncertain  to  report  for  the  higher  sea 
states. 

The  overall  size  of  a  retention  zone  increases  as  the  sea  state 
increases  (since  the  length  of  surface  waves  and  consequently  the  windrow 
separation  and  SRZ  width  increases)  and  as  the  aspect  ratio  increases 
(since  the  cell  penetration  and  consequently  the  SRZ  height  increases). 
Since  the  volume  of  oil  in  a  retention  zone  is  the  oil  concentration 
multiplied  by  the  SRZ  volume,  the  amount  of  oil  trapped  in  an  SRZ 

increases  with  sea  state  and  aspect  ratio.  Notice  that  a  substantial 
amount  of  oil  is  suspended  in  the  retention  zones,  and  these  zones  are 
deep,  even  for  light  winds,  low  seas,  and  relatively  small  LC  penetration 
depths.  Two  points  may  be  relevant  here.  Assuming  the  fractions  of  oil 
shown  in  these  Tables  are  correct,  the  actual  concentrations  of  oil 
(recall  that  our  concentrations  are  normalized  to  give  unit  concentration 
in  the  retention  zone)  in  the  SRZ  will  be  very  small  in  most  cases. 
Secondly,  our  analysis  assumes  that  all  of  the  oil  is  in  the  form  of 

small,  noninteracting  particles.  This  assumption  may  reflect  with 
reasonable  accuracy  situations  in  which  dispersants  have  been  used.  Under 
natural  conditions,  however,  the  oil  to  which  our  work  applies  is  that 
part  which  has  been  torn  off  of  the  floating  surface  oil  layer,  and  one 

must  interpret  the  present  data  as  the  fraction  of  that  oil  is  held  in 

suspension. 

As  a  point  of  reference  in  interpreting  the  results  of  this  section, 
note  that  the  diameter,  d,  of  spherical  oil  droplets  in  water  of 
(molecular)  kinematic  viscosity  v,  will  have  a  terminal  velocity 

Vj  a  1/18  (g/v)(l-po/pw)d2 

if  motion  is  in  the  Stokes  regime.  Thus,  a  one  millimeter  diameter 
droplet  of  oil  with  po/pw  *»  0.99  corresponds  to  Vf  *  0.5  cra/s,  and  if 
po/pw  *  0.98,  Vf  *  1,0  cm/s. 


8.0  Collection  of  Floating  Oil  in  Windrows 

Floating  surface  oil  will  be  gathered  by  Langmuir  circulations  into 
tong  windrows.  This  prominent  feature,  a  generally  reported  one  for  those 
oil  spills  for  which  extensive  observations  are  available  (see  Atwood  et 
al  (1980)),  is  considered  in  this  section. 

The  purpose  of  the  section  is  to  estimate  the  thickness  and  total 
volume  of  oil  in  a  windrow  by  a  method  simple  enough  to  be  used  by  an 
observer  at  an  oil-spill  site.  Oil  is  assumed  to  exist  in  the  surface  as 
a  thin,  coherent  sheet  of  variable  thickness.  The  spreading  pressure  due 
to  interfacial  tension  is  assumed  to  be  zero,  an  assumption  that  is 
reasonable  if  the  oil  is  sufficiently  weathered  and  remains  thick  enough 
to  be  measured  by  mechanical  means.  Under  these  circumstances,  the  oil 
layer  is  influenced  by  gravity,  which  acte  to  spread  the  oil  into  a  layer 
of  uniform  thickness,  and  friction  at  the  oil-water  and  oil-air 
interfaces . 

We  are  interested  here  in  the  build-up  of  oil  held  against  gravity 
by  the  horizontal  sweeping  motions  in  Langmuir  circulations;  variations  in 
the  direction  of  the  wind  are  ignored  (the  stress  imparted  by  the  wind  is 
directed  along  windrows  and  therefore  is  at  right  angles  to  the  direction 
of  interest  here).  Langmuir  circulations  are  assumed  to  exist,  and  to  be 
unaffected  by  the  oil  on  the  surface. 

Adopting  a  conventional  "shallow  water"  model  for  this  problem  (see 
Leibovich  1977c,  Appendix  A,  for  example),  and  assuming  that  a  steady 
balance  between  applied  stress  and  gravity  spread  has  been  achieved,  the 
oil  thickness  varies  with  spanwise  distance  y  (measured  from  the  LC 
convergence  line),  and  rw(y)  the  stress  applied  to  the  oil  layer  by 
Langmuir  circulations  according  to  the  rule 

p0gXh3h/3y  ■  Tw(y)  (8-1) 


where 


X  *  <pw-p0)/P„  (8-2) 

and  g  is  the  acceleration  of  gravity. 

The  shear  stress  on  the  oil  layer  may  be  related  to  the  velocity  (v) 
of  the  water  below  it  by  use  of  a  dimensionless  friction  coefficient, 
cf,  with 
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(8-3) 


Tw  *  (l/2)pwCfV 


Experiments  on  oil  containment  barriers  (Lindenmuth  et  al  1970)  indicate 
that  cf  depends  upon  oil  type  and  is  in  the  range  0.002  to  0.009.  Hoult 
and  Cross  (1971)  adopted  the  (constant)  value  Cf  »  0.008,  and  we  make 
the  same  assumption  here. 

Referring  to  the  sketch  in  figure  (8-1),  suppose  v  is  the  surface 
sweeping  velocity  in  the  LC,  and  LW  is  the  observed  width  of  a  band  of  oil 
collected  in  a  windrow.  Then  h  =  0  at  y  S  LW/2,  and  the  maximum  oil 
thickness,  h^h,,,^  occurs  at  y«  0.  Integrating  (7-1),  we  obtain  the  oil 
thickness  at  any  station  y  as 


h(y) 


f  2  rLW/2 

Lppg  y 


Tu  <*y] 


1/2 


The  total  volume  of  oil  collected  per  unit  length  of  windrow  is 


(8-4) 


Vo  "  2/qW/2  h(y)dy. 


(8-5) 


Suppose  LLC  is  the  distance  between  windrow  centers.  Examination  of 
LC  simulations  shows  that  v  rises  from  zero  at  y  ”  0  to  a  maximum,  say 
vmax>  aC  about  1/3  the  distance  between  convergence  and  upwelling 
planes,  or  about  LLC/6,  and  that  it  falls  back  to  zero  at  the  upwelling 
plane,  which  we  may  put  at  y  “  LLC/2.  Our  estimate  of  oil  collection  into 
windrows  can  only  be  a  crude  one  at  best,  and  it  seems  that  the 
replacement  of  v  by  a  piecewise  linear  form  is  consistant  with  the 
accuracy  of  the  models.  Consequently,  set 


f^Vmax  y/LLC,  0  <  y  <  LLC/6 
v(y)  “  1 

(3v  /2)(l-2y/LLC),  LLC/6  <  y  <  LLC/2 

max  —  —  (8-6) 

Furthermore,  v,,^  will  be  a  fraction  of  the  LC  convective  velocity  scale 
V.  If  we  let  a  represent  the  fraction,  -  aV,  then  our  model  is 


complete.  Inputs  are  LW  and  LLC,  p0t  pw,  a  and  V.  A  survey  of  a 
number  of  LC  simulations  suggests  that  a  *  0.25  is  a  reasonable  estimate 
that  should  be  typical.  The  results  of  the  model  are  being  displayed  in 
dimensionless  form  in  Figures  8-2  and  8-3.  Here  we  have  replaced  V  in 
favor  of  Va  through  (3-18). 

Two  dimensional  cases  are  shown  in  Table  8-1  to  give  the  reader  an 
idea  of  the  numbers  involved;  (3-18)  has  again  been  used  to  replace  V  by 


Table  8-1  -Maximum  oil  thickness  and  volume  of  oil  collected  per 
unit  length  of  windrow  for  a  wind  speed  of  lOms'l,  width  oil  band 
in  windrow  is  assumed  to  be  2  m,  distance  between  windrows  is 
assumed  to  be  5  ra. 


Wind  Speed 

po/pw 

h^axlmm) 

V0(cc/cm) 

5  m/s 

.98 

3.3 

53 

10  m/s 

.98 

6.5 

106 

10  m/ s 

.99 

9.2 

150 

15  m/ s 

.98 

9.8 

160 

As  with  the  other  theoretical  models  in  this  report,  it  is  important 
to  compare  predictions  with  observations  in  the  field;  tuning  of  the  model 
by  appropriate  adjustments  of  constants  is  prudent. 
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Figure  8.2  Maximum  thickness  hp  of  oil  collected  in  windrows.  is  wind  speed;  LLC  is 

distance  separating  windrows:  LW  is  the  width  of  oil  band  in  a  windrow;  g  is 
acceleration  of  gravity;  pq  is  the  oil  density;  pw  is  the  water  density. 


j 


> 


I 


1.0 


V0g1/2 

\yLLC)3/2X1° 


4 


Figure  8-3. 


9. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Our  investigation  indicates  that  oil  is  readily  mixed  in  the  water 
column  to  significant  depths,  provided  the  oil  exists  in  the  form  of  small 
noninteracting  particles.  Oil  is  expected  to  be  in  this  form  when 
dispersants  are  successfully  used;  in  this  case,  the  conclusions  arrived 
at  here  are  encouraging.  Under  natural  circumstances,  a  fraction  of  the 
oil  will  be  torn  into  small  pieces  if  the  seas  are  sufficiently  severe, 
and  this  body  of  oil  will  be  subject  to  the  dispersion  processes  we  have 
analyzed  here.  The  determination  of  the  extent  of  this  prerequisite 
fragmentation  in  specified  environmental  conditions  in  the  ocean  remains 
an  important  open  question. 

In  order  to  make  progress,  we  have  been  forced  to  make  a  number  of 
simplifying  assumptions.  Some  of  these  assumptions  no  doubt  can  be  weak¬ 
ened  or  removed.  The  problem  as  a  whole,  however,  is  of  such  complexity 
that  uncertainit ies  inevitably  will  remain.  It  is  our  feeling  that  the 
present  results,  which  may  not  be  quantitatively  accurate,  incorporate  th  ■ 
essential  physics  and  should  therefore  give  correct  qualitative  trends, 
and  be  correct  so  far  as  orders  of  magnitudes  are  concerned.  Theoretical 
estimates  such  as  ours  can  be  improved  as  data  is  made  available; 
experiments  designed  to  generate  the  required  data  may  be  feasible,  and 
certainly  would  be  desirable. 

The  same  remarks  can  be  made  for  our  simple  model  of  the  collection 
of  floating  oil  into  windrows.  Since  these  oil  bands  are  readily  visible, 
and  the  surface  water  motions  are  monitored  more  easily  than  subsurface 
motions,  it  should  be  relatively  easy  to  "tune"  this  model.  It  has  been 
programmed  on  a  hand  held  calculator,  and  is  therefore  easily  taken  to  the 
field;  if  it  is  verified,  it  should  be  well  suited  for  rapid  estimation  of 
floating  oil  volumes  and  thicknesses. 
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STOMMEL  RETENTION  ZONES 
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Figure  A-4.  See  Figure  4-4 
for  caption. 
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Figure  A-14.  See  Figure  4-4 
for  caption. 


Figure  A-15.  See  Figure  4 
for  caption. 
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Figure  A-18.  See  Figure  4 
for  caption. 
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Figure  A-20.  See  Figure  4-4 
for  caption. 


Figure  A-21.  See  Figure  4-4 
for  caption. 
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